Criticality in Charge-asymmetric Hard-sphere Ionic Fluids 
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Phase separation and criticality are analyzed in z:l charge-asymmetric ionic fluids of equisized 
hard spheres by generalizing the Debye-Hiickel approach combined with ionic association, cluster 
solvation by charged ions, and hard-core interactions, following lines developed by Fisher and Levin 
(1993, 1996) for the 1:1 case (i.e., the restricted primitive model). Explicit analytical calculations for 
2:1 and 3:1 systems account for ionic association into dimers, trimers, and tetramers and subsequent 
multipolar cluster solvation. The reduced critical temperatures, T* (normalized by z), decrease 
with charge asymmetry, while the critical densities increase rapidly with z. The results compare 
favorably with simulations and represent a distinct improvement over all current theories such as 
the MSA, SPB, etc. For z/1, the interphase Galvani (or absolute electrostatic) potential difference, 
A(f>(T), between coexisting liquid and vapor phases is calculated and found to vanish as \T — T c \^ 



when T — > T c — with, since our approximations are classical, (3- 



Above T c , the compressibility 



maxima and so-called fc-inflection loci (which aid the fast and accurate determination of the critical 
parameters) are found to exhibit a strong ^-dependence. 



PACS numbers: 05.70.Fh, 61.20.Qg, 64.60.Fr, 64.70.Fx 



INTRODUCTION 



The location and nature of criticality in ionic fluids 
have been subjects of intense interest in recent years 
HI 13- At sufficiently low temperatures fluid elec- 
trolytes typically undergo separation into low and high 
concentration phases which may be driven primarily by 
the Coulombic interactions. The universality class of the 
associated critical points has been under debate owing 
to apparently conflicting experiments, inconclusive sim- 
ulations, and the analytic intractability of the statistical 
mechanics beyond a mean field level 0, 0, 0] . Possible 
scenarios that have been discussed include, classical or 
van der Waals critical behavior (as might be anticipated 
in view of the long-range Coulomb forces) , crossover from 
classical to Ising-type behavior sufficiently close to the 
critical point [3, Li and, as the leading candidate, three- 
dimensional Ising-type criticality (as might be expected 
for effective short range interactions arising from Debye 
screening): indeed, recent simulations |5j,|g,|2l definitively 
establish Ising behavior for the simplest charge and size- 
symmetric model, namely, the restricted primitive model 
(or RPM); but for z:l and size-nonsymmetric systems, 
the issue is not yet settled. 

The most basic continuum models of ionic fluids are 
the so-called two-component primitive models consist- 
ing of N = N + +N- hard spheres, N + carrying a charge 
q+ = z + qo and N- a charge q- = —Z-qo (with N-/N+ = 
z+/z- = z so that the overall system is electrically neu- 
tral). The background medium is assigned a uniform 



dielectric constant, D, that may be used to represent a 
nonionic solvent. In the simple cases on which we fo- 
cus here, all the spheres have the same diameter, i.e. 

a ++ = a^ =a__ — a. The natural and most appropriate 

reduced temperature variable is then determined by the 
contact energy of a +zq ion with a —qo counter-ion so 
that 



T* — k B TDa^ /q + \q_\ = k B TDa/ zq^ . 



(1.1) 
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Likewise, the normalized density is reasonably taken as 

p* = Na%-/V = pa 3 , (1.2) 

in which V is the total volume. 

This model (with many ionic species) was first ana- 
lyzed by Debye and Hiickel (DH) |£| who derived an ap- 
proximate expression for the Helmholtz free energy by 
solving the linearized Poisson-Boltzmann equation for 
the potential around each hard-core ion. For the sim- 
plest 1:1 (or 2=1) case, i.e., the restricted primitive model 
(RPM), the DH theory predicts [3,L!lj a critical temper- 
ature, T* DH = jq = 0.0625, that is in surpri singly g ood 
agreement with modern simulations |R IfiL FR ITU Il2l ll^| 
that yield T* < 0.05; however, the critical density pre- 
dicted by the DH theory, namely, p* DH = l/647r ~ 0.005, 
is significantly too low since the simulations indicate 
P*c It 0.07. Because ionic criticality occurs at such low 
temperatures, the association of charges of opposite signs 
into 'clusters' is an essential feature in the strongly inter- 
acting regime, as observed in criticality and phase separa- 
tion in other Coulomb systems ^4[- Hence, the first cru- 
cial improvement contributed by Fisher and Levin (FL) 
0, El was to incorporate Bjerrum ion pairing 0] into 
the DH theory: this then depletes the density of the free 
ions that drive the transition, as a result of which the pre- 
dicted critical density increases by a factor of 9. However, 
in order to get an acceptable phase diagram, Fisher and 
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Levin also found it essential to account for the solvation 
of the dipolar ion pairs, or dimers, by the residual ionic 
fluid. The resulting "DHBjDI" theory (with 'DP signi- 
fying dipole-ionic-fluid solvation) yields critical param- 
eters, namely, T* FL ~ 0.055-0.057, p*, FL ~ 0.026-0.028 
which, to date, provide the best agreement with the sim- 
ulations (which indicate T* ~ 0.0493 3 , p* ~ 0.075 @|). 

The other most commonly used theory, the mean 
spherical approximation (MSA) EH E3 EH EH, yields 
T* MSA = 0.0785 (even higher than the simple DH theory) 
and Pc msa = 0.0145 (via the energy route) . Although, like 
the other approximate theories, the FL approach makes 
no reliable statements regarding the universalit y cl ass of 
the criticality — only classical behavior arises [2Cj] — it 
does provide significant physical insights into the origin 
and location of the critical point, specifically identifying 
the role of ionic association, of the solvation of neutral 
clusters, and of excluded-volume effects. 

Two generalizations of the RPM are of profound in- 
terest, namely, the size- asymmetric primitive model (or 
SAPM) and the charge- asymmetric primitive model (or 
CAPM). Indeed, it has been argued f2|, that destroying 
the (artificial) size symmetry of the RPM might even af- 
fect the universality class of the criticality. It may be 
suspected that this feature will eventually be ruled out 
by precise simulations. Nevertheless, it has been demon- 
strated via exactly soluble ionic spherical models [2lj | . 
that size-asymmetry can produce dramatic changes: ex- 
plicitly, the charge correlations become "infected" by the 
critical density fluctuations leading to the destruction of 
normal Debye exponential screening at criticality. Hence, 
asymmetry should be carefully accounted for in any re- 
alistic analyses of ionic criticality. 

In the size-asymmetric model the + and — ions have 
unequal diameters: computer simulations H3] then 
indicate that both T* and p* decrease with increasing 
size asymmetry. However, this is directly opposite to 
the trends predicted by the MSA and some of its exten- 
sions jlil |2J|. On the other hand, a DH-based theory 
developed by Zuckerman, Fisher and Bekiranov j2^| that 
recognizes the crucial existence of "border zones" around 
each ion in which the charge is necessarily unbalanced, 
does, in fact, predict the correct initial trends, as does 
the ionic spherical model [2l| . 

Here we study charge- asymmetric models in which the 
diameters of the basic positive and negative ions remain 
equal, but the charges are in the ratio z:l (z~ = 1). Al- 
though this model is somewhat artificial for applications 
to, for example, multivalent molten salts such as CaF2, or 
AICI3 (since, in actuality, the cation and anion sizes are 
rarely equal), it nonetheless, represents a valuable step 
in searching for a physical understanding of real systems 
[2?j| which exhibit both charge and size asymmetry (as 
well, of course, as other complexities such as short range 
attractions, etc.). 

One should remark, first, that with the normalizations 
tPJ and ljl.21) . the original DH theory Q predicts that 
T*{z) and p*(z) are independent of z; furthermore, the 



same is true for the MSA El El- 

However, we attack the 
problem via an approach which extends the DH-based 
methods developed by Fisher and Levin 0, for the 
RPM as sketched above. Specifically, we calculate ap- 
proximate critical parameters and coexistence curves for 
2:1 and 3:1 systems by explicitly accounting for the asso- 
ciation of the individual ions into dimeric, trimeric, and 
tetrameric neutral and charged clusters, and by includ- 
ing the multipolar cluster solvation free energies induced 
by the ionic medium. In the calculations reported here 
the excluded volume effects associated with the hard-core 
ion- ion repulsion enter in three crucial ways: first, in the 
solvation free energies of the individual ions, as in the 
original DH theory, and of the neutral and charged ionic 
clusters, as in FL; secondly, in the computation of the 
cluster association constants which play a pivotal role; 
finally, via general hard-core 'virial terms' in the free en- 
ergy (described within a simple free- volume approxima- 
tion 

In its primary version our theory may be dubbed 
a DHBjCI approach, with 'CP signifying cluster-ionic- 
fluid solvation including the neutral (z + l)-mer and all 
smaller charged clusters. When specific hard-core (HC) 
excluded-volume virial terms are included, we will la- 
bel the theory DHBjCIHC. More detailed specific re- 
finements will also be examined in order to understand 
the interplay of various effects. However, in all ver- 
sions, our approach unambiguously predicts that the crit- 
ical temperatures, T*{z), decrease with increasing charge 
asymmetry, z, while the critical densities, p*(z), increase 
markedly. This behavior is exhibited in Figs. and 
and clearly contrasts with the z-independence predicted 
by the DH and MSA approximations. Furthermore, one 
sees from the figures that our results mirror closely the 
trends uncovered by simulations [27ll23|. 

The main physical effect behind these trends appears 
to be that increasing the charge asymmetry produces a 
larger number of neutral and charged, but relatively inert 
ion clusters: the depletion of the density of (charged) ions 
and their smaller average mean-square charge leads, first, 
to a lower critical temperature, and, thereby, as in the 1:1 
case, to a higher critical density. In Sec. IIXI we explore 
this interpretation further a nd p resent a comparison with 
other current theories 0, l2fll l3fT | : these either fail to 
yield even the correct sign of the changes with z or else 
predict effects that are much too small! 

In order to obtain efficiently, accurate numerical val- 
ues for the critical parameters implied by the DHBjCI 
theories, we have utilized the so-called ^-inflection loci 
introduced recently [3J, y3|. These are defined as the 
loci on which x — x(T,p)/p k is maximal at fixed T 
above T c , where x{T,p) = p(dp/dp)r- These loci all in- 
tersect at the critical point but their behavior is also of 
interest on a larger scale: See Figs. and 1101 below. In 
our analysis we find that they are strongly dependent 
on the details of the model (such as the hard cores) as 
well as on the charge asymmetry. Thus for our preferred 
parameters, the values of k for which the k- locus has a 
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FIG. 1: Critical temperatures as a function of charge asym- 
metry, z, as predicted by the present DHBjCI theory [see Eqs. 
Hfi.lll . It). 411 and 16. (ill ] and its refinements including hard-core 
(HC) virial terms with 'standard' (triangles: see Table |n|l and 
'optimal fit' parameters [crosses: see Eqs. I|fi.2l) . lfi.5| and 
lfi.9l ]. compared with Monte Carlo simulations [27ll2g|| (open 
circles) and the original Debye-Hiickel (DH) theory. The spe- 
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when z ^ 1, This electrostatic potential, appropriately 
deemed a Galvani potential l34| . has been explicitly 
anticipated in the case of 1:1 electrolytes with nonsym- 
metric ion-ion interactions in an interesting phenomeno- 
logical treatment by Muratov j^| (that, however, fails to 
satisfy the important Stillinger-Lovett sum rule [3J). It 
also features in a detailed discussion of colloidal systems 
(with z > 1) by Warren |37j . 

However, the dependence of Atfi on z for moderate 
charge asymmetry has not been examined previously. On 
approach to the critical point we find that A0(T) van- 
ishes as \T-Tcl 13 , where, because our approximations are 
classical in character, we obtain /3 = s. (A similar conclu- 
sion is reached for the asymmetric 1:1 electrolyte in 35]). 
As a consequence of this potential difference, a charged 
double layer [sH. 03| will exist at a two-phase, liquid- 
vapor interface; this, in turn, will be of significance for 
interfacial properties such as the surface tension, which 
have been studied theoretically, to our knowledge, only 
for the RPM jjjUjjjJ. 

The balance of this article is laid out as follows: in 
the next section pertinent thermodynamic principles are 
summarized briefly. Sec. lIIII then describes the computa- 
tion of association constants for the "primary" set of as- 
sociated clusters consisting of one cation of charge +zq 
and m < z anions of charge —go; detailed calculations 
for tetramers are presented in Appendix El The crucial 
multipolar electrostatic contributions to the Helmholtz 
free energy are analyzed in Sec. IIVI These and other in- 
gredients are combined in Sec. to obtain expressions 
for the total free energy and, thence, in Sec. IVII quan- 
titative results for the 2:1 and 3:1 models (following a 
brief account of the pure DH theory). A discussion of 
the ^-inflection loci is presented in Sec. IVIII Sec. IVIIII 
is devoted to the Galvani potentials while, as mentioned, 
our results are reviewed briefly and compared with those 
of other current theories in Sec. lIXI the varied predictions 
for the critical parameters being summarized in Figs. 1131 
andGJ 



II. SOME BASIC THERMODYNAMICS 



A. Phase Equilibrium 



FIG. 2: Critical densities as a function of charge asymmetry 
z, as predicted by the DHBjCI theory and its refinements, 
using the same conventions as in Fig. ^ 



vertical slope at (T c ,p c ) are k (z)~ 0.93, 0.18, and —0.87 
for 2 = 1, 2, and 3, respectively. It should be possible to 
check these predicted trends via simulations. 

An interesting new feature that arises in our calcula- 
tions (but is absent by symmetry in the RPM) is the 
appearance of a nontrivial electrostatic potential differ- 
ence, A</>(T), between coexisting liquid and vapor phases 



A z:l electrolyte may be regarded (neglecting the sol- 
vent) as a single-component system since putting Nq 
neutral 'molecules' (each of one positive and z negative 
ions) at temperature T into a domain of volume V com- 
pletely defines the thermodynamic state. The total num- 
ber of ions is then N = (z + 1)Nq, while the total ionic 
number density, p = N/V, also measures the density of 
the original molecules. The total Helmholtz free energy, 
F(T, V, N) may be introduced, in standard notation, via 
the differential relation 

dF = -SdT -pdV + fidN , (2.1) 
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where p is the chemical potential conjugate to the total 
number of ions. In the thermodynamic limit, the reduced 
variables 

f{T,p) = -F/Vk B T, and 

p(T, p) = p/k B T = -{df/dp) T , (2.2) 

are convenient 0,^3 . The reduced pressure follows from 
the variational expression 

p(T, n) = p/k B T = max (/(T, p) + pp) . (2.3) 

p 

Then phase coexistence (if present) at a given temper- 
ature is specified by the equilibrium conditions 

P(T, pv) = p(T, pi) and p(T, p v ) = p(T, pi) , (2.4) 

where the subscripts v and I indicate vapor and liquid 
phases, respectively. These equations determine the den- 
sities in the two phases: at the critical temperature and 
density, pi{T) and p v {T) coincide. 

The single-component or molecular thermodynamic 
formulation takes care of overall electroneutrality in a 
natural way and utilizes only one overall chemical po- 
tential. It is complete, in principle, if one knows the 
Helmholtz free energy density f(T, p) . An alternate ap- 
proach is to treat the "isolated" or free ions, and the var- 
ious clusters into which they associate, each as distinct 
species in thermal equilibrium with one other. Since the 
exact calculation of f(T,p) is intractable, this latter ap- 
proach is useful in constructing approximations for the 
overall free energy density. Such a formulation, however, 
requires the principles of mw^icomponent thermodynam- 
ics that have have been reviewed s yste matically by FL 
for the charge-symmetric RPM in (henceforth ab- 
breviated as I). The formulation needed for the charge- 
asymmetric models is quite similar to that outlined in I 
but contains some subtle differences. Thus, even at the 
cost of some repetition, we outline the main principles 
here. 

Consider a system of distinct species er, which may 
be free ions or ion clusters (a = +, — for the original 
ions, and a = 2, 3, • • • for dimers, trimers, etc.), with 
number densities p a = N a /V, where is the number of 
entities of species a. The Helmholtz free energy density 
/(T; {pa-}) can be defined through a generalization of the 
single-component formulation above [see Eqs. (2.4) and 
(2.5) in I]. The reduced chemical potential for species a 
then follows from 

p a (T; {p a }) = p a /k B T = -(df/8p a ) . (2.5) 

Since all the species present will be in chemical equilib- 
rium, the sum of the chemical potentials of the reactants 
in any reaction will equal the sum of the chemical po- 
tentials of the products [see 1(2.2) and 1(2.3)]. These 
equations together with conditions l|2.4l) and overall elec- 
troneutrality, namely, 

]T q aP „=0, (2.6) 



determine the system in equilibrium. For calculating the 
pressure one may still use Eq. H2.3I) . or, equivalently, the 
multicomponent form 1(2.6). 

For a multicomponent system in which none of the 
species has a net charge, thermal equilibrium demands 
that the chemical potentials of each species match in co- 
existing phases. More generally, however, it is the elec- 
trochemical potential that must be equal in both phases 
so that for a species a one has 

P<r,v + q<r4>v = P-a,l + q<j<t>i , (2.7) 

where q a is the net charge of particles of species a and 
4> V (T) and 4>i{T) denote the (in general distinct) elec- 
trostatic potentials in the coexisting vapor and liquid 
phases, respectively. Then A<j)(T) = <pi — <p v is the ab- 
solute electrostatic potential difference between the two 
phases, i.e., the interphase Galvani potential (33l l34l|. In 
the molecular or 'overall' formalism presented above, the 
correct phase behavior can be obtained without any ref- 
erence to the Galvani potential since the chemical poten- 
tial p, conjugate to the overall density p, corresponds to 
a neutral species that is insensitive to the electrostatic 
potential <\>. Nevertheless, the Galvani potential repre- 
sents a significant feature, that is not present (or vanishes 
identically) in the RPM: it is discussed in further detail 
in Sec IvTtTI 



B. Free Energy Contributions 

Our aim is to construct a physically appropriate, al- 
beit approximate free energy for the model systems by 
adding contributions that arise from the various degrees 
of freedom and the underlying mechanisms and interac- 
tions. As a zeroth order approximation any fluid may be 
taken as an ideal gas. Thus, for each species we invoke 
an ideal-gas term 

f d (T; p er ) = Pa - Pa HfMT)) , (2.8) 

where C a (T) depends on the internal configurational par- 
tition function of species a, C<?(T), and the de Broglie 
wavelength, A CT (T) (see I for details). 

The principal contribution to the interaction free en- 
ergy of our model electrolyte comes from the electrostatic 
interactions between the ions. We will use a DH "charg- 
ing" approach to calculate the electrostatic free energy of 
each species as discussed in detail in Sec. IIVI below. The 
only other significant interaction between the ions is the 
hard-core interaction. 

The various forms of additive free energy corrections 
for the hard-core contributions that might be employed 
are discussed at length in I. However, we have not ex- 
plored the range of these options here. It may be noted, 
first, that such second-order and higher virial-type cor- 
rections [9, LLJj] ; enter formally in higher order in powers 
of the overall density than do the (leading) electrostatic 
terms. Secondly, the exact hard-core diameters already 



5 



play a quantitatively significant role in DH theory itself 
(see Sec. IV Al below V Furthermore, as observed in the 
Introduction, the exact hard-cores are equally vital in 
the formation of ion clusters, thereby affecting the values 
of the corresponding association constants which in turn 
play a dominant role. Finally, the formation of tightly 
bound clusters (see I and below) at temperatures < T c 
has the effect, at the rather low densities near criticality, 
of markedly increasing the available free volume relative 
to a fluid with only hard-core interactions. 

For these reasons, in the present study we have con- 
fined our considerations to a simple free-volume approx- 
imation which adds 

/ HC = (E /OMl-£ B <rP«), ( 2 -9) 

to the Helmholtz free energy. In the low density limit 
with all species regarded as hard spheres of diameters 
a CT , the exact value of the coefficients B a is 27ra^./3. But, 
as noted in I, this choice for equisized hard spheres im- 
plies an unrealistically low maximum density at p* — 
3/2ir ~ 0.48 in contrast to the true, fee packing density 
°f Pmax = V 7 ^ ~ 1.41. A reasonable alternative choice 
for use at intermediate densities is thus to take effective 
values of the B a coefficients corresponding to bec close 
packing H}, namely, B a /al = 4/3^3 ~ 0.770. We will, 
hereafter, refer to this choice as using "bec hard cores"; 
its influence on the values of the critical parameters will 
be examined below in Sec. IVII 

We may remark, however, that while some improve- 
ments in accounting for volume-exclusion may still be 
feasible, the studies in I indicate that straightforward, 
naive approaches tend to strongly over-estimate the ex- 
cluded volume effects. This seems to occur because the 
dominant many-particle ionic correlations in the low- 
temperature moderately dense liquid lead to an "ex- 
panded crystal-like" structure that screens out direct 
hard-core interactions: see plot (d) in Fig. 6 of I and 
the related discussion in I Sec. 8.5. 



III. ASSOCIATION CONSTANTS FOR ION 
CLUSTERS 

The success of the DHBjDI theory in estimating the 
critical point of the RPM, together with "snapshots" of 
primitive models of ionic fluids from computer simula- 
tions [27l l2l| indicate that a high degree of association is 
present in the critical neighborhood. A careful analysis of 
the association constants for ion clusters is therefore es- 
sential. Here, we generalize Bjerrum's original approach 
[lfl| to define and calculate the association constants for a 
set of "primary" clusters which contain one central cation 
of charge zqo surrounded by 1 < m < z singly charged 
anions (a general dimer configuration is illustrated in 
Fig. El. These primary clusters with, including the bare 
ions, net charges q a = (—1, 0, +1, ■ • • , +z)q n , will be the 
first to form when the temperature is lowered, and the 
density, increased. 



z 



i 




FIG. 3: A dimer with the two ions separated by a distance ai . 
The dotted sphere indicates the closest possible approach by 
a screening ion. The dashed sphere of radius a 2 represents an 
effective exclusion zone for solvation computations: Sec. lIV Al 



Of course other, larger clusters of ions must eventually 
come into play. However, insofar as their net charges 
fall in the same range, they may be subsumed, in a first 
approximation, under the like-charged ions and primary 
clusters: see the discussion in I Sees. 8.3 and 8.4. The 
most important exception is probably the doubly over- 
charged "molecular" clusters with m — z + 2 counterions: 
however, we believe that such clusters will not contribute 
significantly in the critical region for z < 4. 

To proceed, consider a charge q+ = zqo fixed at the 
origin of a Cartesian coordinate system with m satellite 
charges q-=—qo around it. For m = 3 one has a tetramer 
for which Fig.Q]illustrates a general configuration, and let 
Ti be the position vector of the i~ th satellite. The reduced 
configurational energy (electrostatic plus hard-core) for 
such a system, normalized by q + q-/Da, is 

E m ,z({ v i\) = — X! — ' if r *» - a > 

i=l (1,2) 

— —00, otherwise, (3.1) 

where r, = |r,|, r.^ = |r^ — rj\, and (i, j) indicates a sum 
over all distinct pairs. The association constant for a 
cluster or (m + l)-mer formed by these charges is the 
internal partition function 

K m , z (T; R) = —J[ dri exp[E m , z ({Ti})/T*] , 
m ' i=i-' a 

(3.2) 

where R is a suitable cut-off radius without which all the 
integrals would diverge at large distances. 

The choice of R is necessarily somewhat arbitrary since 
there is no clear, absolute criterion for when a group of 
ions is to be considered "associated". The ambiguity in 
choosing a cut-off radius arises even for the simplest pos- 
sible cluster, a dimer (m—1). In that case Bjerrum ^| 
observed that the integrand in 113.211 exhibits a minimum 
at a radius R B] = a/2T*; thus he chose R = i? Bj as the 
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FIG. 4: A configuration of a tetramer with coordinates suit- 
able for calculating the association constant. The small dot- 
ted spheres indicate the ground-state orientations of the satel- 
lite ions. The central positive ion is located at the origin. 



cut-off for T* < \ (since for T* > § one has R B > < a). 
Evidently, the choice R = R Bl , makes the association con- 
stant least sensitive to the value of the cut-off. Bjerrum's 
choice, however, may reasonably be considered as un- 
physical since R Bi (T) becomes unbounded when T falls 
to zero while one expects a (+, — ) ion pair to become 
more tightly bound at lower temperatures. This issue 
is discussed in I (see Sec. 6.2) where Bjerrum's associa- 
tion constant is also compared with other definitions: see 
also [25L l4lj | . While Bjerrum's cut-off, R B ', has no direct 
relevance to the actual physical size of a dimer which is 
much more compact (as analyzed in I, Sec. 6.2) — the 
value and behavior of Bjerrum's association constant is 
numerically accurate for T* < despite the unphysical 
nature of the cut-off. 

In light of this analysis we generalize Bjerrum's ap- 
proach and choose the cut-off R so that (dK„ uz /dR) is 
minimal. For dimers, the choice of the cut-off makes 
very little difference over a wide range of R at and below 
T* = j£. However, one may anticipate that the depen- 
dence of K m=ZtZ (T; R) for to> 1 will be more sensitive to 
the choice of R because the ground state binding energy 
per (5+, q-) bond (for a neutral cluster) becomes smaller 
with increasing z ■ As a consequence we must expect 
our estimates for the critical densities to become less re- 
liable with increasing charge asymmetry. 

It should also be noted that our choice of integration 
domain H3.2JI is somewhat arbitrary. Thus by taking the 
outer boundary surface to be n = R (for all i) we have 
chosen to integrate over an m-dimensional hypercube. 
Instead, one might well choose the m-dimensional hyper- 
sphere: J2j r j < R 2 + (m — l) 2 a 2 (where j runs from 1 



to to), or, say, a hypercube cut along its body diagonal, 
J2j r j < R + (m — l)a. However, it is reassuring that 
for to = 2, where exact numerical calculations are possi- 
ble, the choice of integration domain makes a difference 
of less than 0.5% in ^2,2 at T C * DH = 0.0625; furthermore, 
the sensitivity to this choice is reduced at lower temper- 
atures. 

Now, for small T, the integral defining K m ^ z is domi- 
nated by the ground state energy of the cluster, and, it 
is appropriate, therefore, to expand the integrand about 
the ground state configuration. After appropriate scaling 
of the radial variables, we obtain, for to > 3, the general 
form 



K m . z (T;R) 



1 8ir rn+1 / 2 J m a 



3 m 



to! rr 2 ™- 3 \V 2 

z m-3/2 /jn*2m-3/2 



1JU=1 A m,k 

eMmC m , z /T*)X mtZ (T*;R), (3.3) 



where the residual integral satisfies 

l m , z (T*;R) = l + 0(T* 



(3.4) 



while the Xm.kS are the eigenvalues of the reduced 
quadratic form describing the angular variation of the 
energy. The dominant exponential dependence is con- 
trolled by C mtZ , the binding energy per satellite in units 
of (q+q-/Da), while J m is the Jacobian of the trans- 
formation leading to the diagonalization of the angular 
integrals. In Appendix El the calculations are performed 
explicitly for tetramers (to = 3: see Fig. thereby il- 
lustrating the general procedure. Evidently the principal 
T-dependence of the association constant can be found 
by scaling the variables and expanding the integrand for 
small T*. A full calculation, however, requires an evalu- 
ation of the residual integral factor X m . z (T*;R). 

Dimers and trimers turn out to be special cases for 
which the general form l|3.3l) does not apply. Neverthe- 
less, the calculations follow similar lines. For dimers the 
association constant has been discussed in detail in I. 
Expanding around the ground state yields 

K lt ,(T; R) = 4™ 3 T* exp[l/T*]I llZ (T*; R) . (3.5) 

Moreover, using the Bjerrum cut-off R Bi and evaluating 
analytically the integral over r in ll3~2l gives [13] 



l hz (T*;R B >) - 



6T* 



r 1/T * [Ei(l/T*)-Ei(2) + e 2 ] 

' (1 + T* + 2T* 2 ) , (3.6) 



6T 



where Ei(y) is the standard exponential integral. Because 
of the normalization (1.1), this result is independent of 
z. The asymptotic expansion for small T* is, in addition, 
independent of R and given by 

l hz (T*;R) = l + 4T* + 4-5T* 2 +4-5-6T* 3 + ... . (3.7) 



7 



We note that when T < 0.1, this expansion gives reason- 
ably accurate results if truncated at the smallest term: 
see I. 

For trimers, a similar but more elaborate calculation 
yields, 



K 2>Z (T*;R) = 32tt 



zT* 3 a 6 



(1 - l/4z) 2 
exp[2(l-l/4z)/T*]r 2 ,,(T*;i2). 



(3.8) 



However, in contrast to dimers, an exact analytical re- 
sult for 22,* (T*; R) seems inaccessible for any value of z. 
Nevertheless, one can obtain precise results by numerical 
integration. For our subsequent calculations we need re- 
sults for trimers with z = 2 and 3, i.e., the integrals 2 2j2 , 
and 22,3- Generalizing the Bjerrum procedure, we deter- 
mine the appropriate, optimal cutoffs, R m ,z by search- 
ing numerically for the minima of (dK m ^ z /dR) at fixed 
temperature. The results can be scaled conveniently by 
setting R2,z/a — R2,z{T*) /T* , where, typically, we find 



7? 2j2 (0.05) ~ 0.263, ^,3(0.05) ~ 0.336 . 



(3.9) 



For z = 2, the sensitivity of the trimer association con- 
stant to the cut-off is markedly greater than found for 
dimers (which was illustrated graphically in I Fig. 3). At 
a temperature T* = 0.052 (some 6% above the predicted 
value of T*(z = 2): see Sec. ED increasing (decreasing) 
Ri,2 by 20% increases (decreases) K2.2 by about 0.7%. 
However, as is explained in Sec. IVII these changes do 
not significantly alter the predicted critical parameters. 
Because of the larger central charge, the sensitivity of 
7 2 ,3 to the cut-off is significantly smaller at the relevant 
temperatures. 

Even if one employs a precise numerical calculation of 
K2, z for trimers, it is worth noting that it can be repro- 
duced quite accurately by Pade approximants that 
embody the small-T* expansion of 2 2jZ . For example, 
when 2 = 2, the expansion 



X 2 , 2 (T*;i?) = 1 



76 
49 



222368768 3 
117649 



357248 2 

2401 
7109382144 4 

117649 



(3.10) 



yields a [1/3] Pade approximant that up to T* = 0.055, 
agrees with the results of numerical integration to better 
than 4%. With this in mind, we consider the approach 
also for tetramers. 

Indeed, for tetramers (needed only for the case z = 3), 
expanding about the ground state yields, 

K 3 , Z (T*;R) - ffiw/V/' f /2 f >T *° /2 3 x 

5 (1-1/V3z) 3 
exp[(3- V3/z)/T*]l 3 , z (T*;R); (3.11) 

See AppendixEl However, the calculation of 2 3-z proves 
difficult even numerically. Asymptotic expansion for 



TABLE I: Fitted expansion coefficients ij for calculating the 
tetramer association constant K 3 , 3 (T*): see ft-i.llfl - 13.1311 . 



3 




3 


10- (j+3) ij 


3 


lO" 16 ^ 


8 


- 0.419627 


12 


55.247 


16 


13.8829 


9 


2.17887 


13 


- 44.2769 


17 


2.58295 


10 


2.6276 


14 


12.4183 


18 


0.40201 


11 


- 28.3178 


15 


0.558599 







small T* 
efforts) 

r(7) 



yields 2 



3.3 



_ T (7) 
— x 3,3 



0(T* 8 ) with (after some 



X^(T*;R) = 1 + 4.26324 T* + 157.697 T* z 
+ 353.407 T* 3 + 29636.117T* 4 - 58642.1 T* 
+ 8.5259.10 6 T* 6 - 7.07815. 10 7 T* 7 . 



(3.12) 



One can then form and examine all the approximants up 
to order 7. One observes readily that the [5/2] approx- 
imant seems the most reliable judging its convergence 
relative to the other approximants: see Fig. 

However, since the tetramer is of prime importance 
for criticality in the 3:1 model, and because one knows 
that an approximant based only on the low-2 asymp- 
totics must fail at some value of T* (of likely magnitude 
~0.11we have undertaken a Monte-Carlo evaluation of 
23,3 0| : see Fig- The details are described in Ap- 
pendix El It transpires that the [5/2] Pade approximant 
agrees to within 4% with the precise numerical calcu- 
lation up to T* = 0.06 (which is 15% higher than the 
DHBjCI value of T*(z = 3) as can be seen in Fig. 
Nevertheless, for our explicit calculations we fitted the 
Monte Carlo calculations of 23 3 to the form 



2 3 ,3(T*; J R)=2< 7 3 ) (T*;i?) 



V 18 ij T*i 

^j=8 J 



(3.13) 



where the coefficients ij are listed in Table ||| As seen 
in Fig. the fit is very good and, indeed, provides an 
accuracy of one part in 10 3 or better. In reality, it prob- 
ably remains valid some way above T* = 0.10; but it is 
also clear that the [5/2] approximant fails rapidly above 
T* = 0.06 and shows noticeable deviations already for 
T*>0.04. 



IV. ELECTROSTATIC CONTRIBUTIONS TO 
THE FREE ENERGY 

A. General considerations 

To calculate the electrostatic part of the free energy we 
adopt the basic DH strategy y| but, as in Q, we gener- 
alize the approach to include cluster species that contain 
more than one ion and, thus, are not spherically sym- 
metric. Consider a cluster (possibly just a single ion) of 
species a that has charges {qi} at positions {n}. Ow- 
ing to the hard-core repulsions the "free" screening ions 
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FIG. 5: Calculation of the association constant integral 
1z, 3 (T*). The dashed line represents the [5/2] Pade approx- 
imant while the dotted lines portray the [6/1] and [2/5] ap- 
proximants. The solid circles result from Monte Carlo inte- 
gration while the solid line is a polynomial fit: see l|3.13fl . 



are prevented from entering the "exclusion zone" of the 
species: see, for example, the dimer with one +2qo and 
one — <?o ion which has a dumbbell shaped exclusion zone 
as seen in see Fig. 

To estimate the free energy of an isolated cluster in an 
atmosphere of screening ions, of densities p v and charges 
q v , we approximate its exclusion zone by a sphere of ra- 
dius do- jjj, El : for the selection of an appropriate value 
for a a , see below in Sec. IV CI At this point we will sup- 
pose only that the choice of origin for this effective ex- 
clusion sphere is such that all the charges of the cluster 
are included within it. A specific criterion for the precise 
choice of origin for the effective exclusion sphere (when 
not dictated by an obvious symmetry) will be developed 
for each cluster species as we address them individually. 

For r < a a , the overall electrostatic potential may gen- 
erally be expanded in terms of the spherical harmonics 
Y\ m as 



4tt 



D ^ 21 

l,m 



^ 9i*Lm(0;><ft)- v 



j+i 



A, 



Y l>m (9,<p), (4.1) 



where D is the dielectric constant of the medium, i 
labels the particles of the cluster a, the are their 
charges and the (ri,6i,tpi), their coordinates, while r^< = 
min(r, r^), r, ( > = max(r, r,), and the notation J2i m 

means S^oEm=-i- We note that the boundary con- 
dition at the origin is already taken into account in this 
expression. 

For r > a a the potential arising from the cluster ions 
in a is screened by the external ions and hence we may 



expand the potential as 

*>(r,0,<p) = ^X>, m fci(Kr) Y l>m (6,tp), (4.2) 



in which screening is embodied in the boundary condition 
<J>> — > when r — > oo [which relates rather directly to the 
introduction of the electrostatic potential 4> in 12.71 ]. The 
inverse Debye length introduced here is defined generally 
by 



k(T,{ Pt }) = (47rJ2 T pTq 2 T /Dk B T 

and, when convenient, we will write 

na = x and na a = x 
The spherical Bessel functions 

h(x) =gi{x)e- x /x l+1 



1/2 



1/£d, (4.3) 



(4.4) 



(4.5) 



that arise in the solution of the Debye-Hiickel or lin- 
earized Poisson-Boltzmann equation, are conveniently 
specified in terms of the polynomials 



9l 



(*) = E 



(l + m)\ ^ 



2 rn m\(l - to)! 

m— v 7 



(4.6) 



(so that go{x) = l, gi(x) = l + x, g2(x) = 3 + 3x + x 2 , etc.) 

On the surface of the exclusion sphere, r = a a , match- 
ing $< and V<i>< to and V$> (the usual conditions 
expressing continuity of the potential and absence of sur- 
face charge) yields the coefficients 



A 



l.m 



Q. 



l.m 



21 + 1 



B, 



(2l + l)h(x a ) 
x a h +1 {x G ) 



4~nQl, m 
a^ 1 X a h+l(x a ) 



(4.7) 



(4i 



in which the cluster multipole moments, Qi, m , which will 
play a central role in our calculations, are given by 



Qi, m = y2.Y* (9i, ^ q t n l , 



(4.9) 



where the summation runs over the particles of the clus- 
ter er. In <4.H . the potential arising directly from the 
ions in the cluster (without any contribution from the 
screening ions) is 

X>l£ m (0i,¥>i)%T> (4.10) 



l.m 



9 



and therefore, the potential inside the exclusion sphere 
arising from the external screening ions is merely 

*<(r;Kft}) = E2|TiAm^l,m(^). (4-11) 

l,m 

where r = (r, 9, ip) and the A^ m are given by 114.711 . The 
electrostatic contribution of the species a to the total free 
energy now follows via the Debye charging process 0, 
as 

AT /"l 

(4.12) 

and normalizing by Vk s T, we finally obtain 



(3 ^ 



47T 1>2!(X CT ) 



|2 

m I 3 



(4.13) 



where the crucial expressions are 

(2^ + l)fc ; (Ax) 
Ax hi+i(\x) 



V2i(x) = d\\ 
Jo 

21 + 1 



1 - 



In 



5/+i(0) 



2(2/ + 1) 



(4.14) 



while the 'multipole-squared amplitudes', J2 m \Qf.m\ 2 
are independent of the axes defining the polar coordi- 
nates. 



B. Monomers 

Consider a monomer (or single + or — ion) with di- 
ameter a±=a and charge q± . The multipole expansion 
(14, 9 j) contains only the I = term with Qo,o = 9±/V47r. 
Substituting into H4.13I) gives the reduced free energy of 
a monomer in a cloud of screening ions 



fl\T,{ Pa }) = 



(4.15) 



If only monomers are present, summing the contributions 
from the positive and negative ions leads to the familiar 
DH free energy 0,0, namely, 

f DH (T, {p a }) = [ln(l + m) - na + i( K a) 2 ]/47ra 3 . (4.16) 

This result, which depends only on x = na is, in fact, gen- 
erally valid for any number of charged species, provided 
all of them have the same size and the system is over- 
all electrically neutral: as well known, it reproduces the 
exactly known answers at the leading low-density order. 



C. Dimers 

For our z:l system, consider the dimer illustrated in 
Fig. with a cation of charge q + = zq , separated from 
an anion of charge g_ = — go by a distance a\. (In reality, 
ai > a will be a fluctuating distance; but, as discussed in 
detail in I, we may, in reasonable approximation, regard 
it as a definite function of T: see also in Sec. IV CI be- 
low.) Let et2 be the radius of the effective exclusion sphere 
(which, clearly, should increase when a\ increases). Since 
the dimer is asymmetric unless z = 1 we displace the cen- 
ter of the exclusion sphere towards the positive cation, 
by a distance pay. see Fig. 03 One should expect the 
optimal value of p to depend on z: by symmetry, one 
must surely choose p=\ f° r z = l> but when z — > oo one 
should, likewise, have p — > 0. For the moment it suffices 
to assume < p < 1 : a concrete criterion for choosing p 
will emerge below. 

In the configuration of Fig. the leading multipole 
moments are 

Ql,0 = \l^-q a[[zp l + (-ly+^l-p) 1 ] , (4.17) 

and Qi t m — if m ^ 0. Substituting the above into 114. 13|) 
yields the dimer contribution 

ir(T;{p ff }) = ^x 

oo 21 

The value of this sum and its rate of convergence clearly 
depends on the value of p. Explicit numerical tests using 
a\/a = l, aija = 3[1 + l n(3)/ 2]/4 (the 'angular average' 
value discussed in Sec. IV Cll show that the series con- 
verges sufficiently rapidly that, to the precision of inter- 
est, one need not consider terms beyond 1 — 2 (see also I). 
Indeed, for z — 2 and 3 and 1 < X2 < 5, the l>3 remain- 
der for p= i varies only from 0.8% to 1.6% of the I = 2 
or dipolar term and can thus be safely neglected within 
the accuracy of this calculation. Evidently, a reasonable 
criterion for the optimal value of p would be that which 
minimizes the full sum of terms from I = 2 to oo. In view 
of the rapid convergence, however, a very satisfactory op- 
tion is to choose the value of p that minimizes the I — 2 
or quadrupolar term: this yields the simple result 

(4.19) 



p = l/(l + Vi) 



This value, in fact, eliminates the quadrupolar term en- 
tirely and satisfies the two limiting cases, z = 1 and 
z — > oo, discussed above. Adopting this expression for 
p and neglecting the terms with I > 3 in l|4.18jl . we ob- 
tain the very satisfactory approximation 



mT ] { Pa }) = 



(z-l) 



T , P2—Vo(x 2 ) 

zl * a 2 



T 



1 aai , 

P2^-V 2 (X2), (4.20) 
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FIG. 6: A trimer bent at an angle 2a. The dotted sphere 
indicates the closest approach by a screening ion. The ex- 
clusion zone is approximated by a sphere, shown dashed, of 
radius 03. 

which we will employ below. It is interesting to note that 

the choice Q4.19J1 also makes the coefficient of the I = 1 (or 

dipolar) term independent of z. Furthermore, numerical 

calculations show that for this choice of p, the I > 3 

remainder divided by the I = 1 term is reduced by a factor 

of about 1 lz relative to the symmetrical assignment p = 
1 

2 • 

D. Trimers 

In considering the solvation of a trimer species, the 
first point to note is that although the ground state is 
linear (in the form: — qo, +zqo, —qo) and so has a vanish- 
ing dipole moment, the typical fluctuating configuration 
at finite temperatures must be bent and hence have a 



dipole moment of magnitude of order q$a. Indeed, ex- 
amination of snapshots of simulations for z > 2 in the 
critical vicinity (see, e.g., [13]) fully confirms this con- 
clusion. Accordingly, consider, as illustrated in Fig. a 
trimer which is (say, "instantaneously") bent at an angle 
2a. To simplify the analysis, we will suppose that it is 
adequate to fix the radial distances r\ and r 2 for both 
satellite anions at the spacing 01. (As discussed in I, 
and also below, we expect the fluctuations in 01 to be 
relatively small.) 

For the effective exclusion sphere, now of radius, say, (13 
(see Fig. 01 , the issue of the placement of its center again 
arises. By symmetry (having imposed ?'i = r 2 =ai) the 
center should lie on the bisector of the angle 2a which, in 
Fig-El nas been identified as the z-axis. Then, in analogy 
to the dimer, we center the exclusion sphere at a distance 
paicosa displaced from the center of the central cation 
(or charge q+) towards the two anions of charge q- whose 
axial location lies at a distance a\ cos a as projected onto 
the bisecting axis. With this placement of the center, we 
find the multipole-squared amplitudes 

|Q ,o| 2 = (l/4^(z-2) 2 g 2 , (4.21) 
]T \Qx, m \ 2 = (3/4tt)[*p + 2(1 - p)} 2 (cos a) 2 q 2 al, 

* '771 

(4.22) 

V \Q2, m \ 2 = (5/47r){3sin 4 a+ (sin 2 a 

+ [^ 2 - 2(1 -p) 2 ] cos 2 c^a 4 . 

(4.23) 

In an ideal calculation of the solvation free energy of 
trimers, every trimer bent at a specific angle would be 
treated as a separate species in its own right. However, 
to make our calculations tractable we substitute these 
expressions for the multipole-squared moments into the 
basic result 114. 13|) and replace the factors that depend on 
a by thermal averages to obtain 



fi\T;{p a }) 



P3 

zT* 



(z-2) 2 -vo(x 3 )+lzp + 2(l~p)} 2 < 



'COS Q 



«3 



/ -3-^2(^3) 



4 



+ (3 sin 4 a + {[zp 2 ~ 2(1 - p) 2 } cos 2 a + sin 2 a} 2 ) — f v A (x z ) + 



(4.24) 



Again, an ideal calculation would recognize that the in- 
creased solvation free energies resulting from larger dipole 
moments, should enhance the thermal weight of more 
highly bent trimers. However, we will forgo such a refine- 
ment (which would require a cumbersome self-consistent 
formulation) and merely weight the bent trimer config- 
urations via the Boltzmann factors computed with the 
"bare" cluster energies. Accordingly, we consider the 



thermal average (O) of an angular function O at tem- 
perature T to be defined by 

(0)= da sin2aO(a)e £(Q)/T 7 

I da sm2ae E( - a)/T ' , (4.25) 

Jtt/6 
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where E(a) = — a/2za\ sin a is the reduced repulsion en- 
ergy between the two satellite anions. Note that in set- 
ting the lower limits of integration at a = ir/6 (due to 
hard-core repulsions) , we have neglected a domain of clos- 
est approach, and, hence, highest repulsive energy, that 
is accessible when ai > a. In fact, in the following calcu- 
lations, we will need only the two averages (sin 2 a) and 
(sin 4 a) which follow from the expressions 

an v = 1 s 2n+3 (l/zT*)- S2n+3 (l/2zT*) 
[ a> {2zT*Y n s 3 (l/zT*)-s 3 (l/2zT*) 

(4.26) 

where 

s n (x) = \ L \.. Ei(-x) + — V(-l) fc ^ . 
w (n-1) x ^ x k 

(4.27) 

As regards the choice of p, a first guess is to choose the 
value Pmin which minimizes the quadrupolar term. How- 
ever, this leads to unphysical features such as p m m < 
and even to p m in — ► — oo when T* — > (when, in fact, 
trimers become straighter and straighter). The alterna- 
tive adopted here is to accept the value of p in the inter- 
val < p < 1 that minimizes the quadrupolar term. One 
may verify that this value is p — for 0.003 < T* < 0.06, 
yielding a quadrupolar term that agrees with the exact 
minimum to within 3%. 

In summary, although, as indicated, various refine- 
ments of our approach may readily be contemplated, 
we believe that the formulation reasonably captures the 
essential physics underlying the solvation of fluctuating 
trimeric ion clusters. 



E. Tetramers 

For z > 3, one must allow for the formation of 
tetramers and include their solvation free energy. A 
tetramer in its ground state is planar with <p=0, satellite 
radii r; = jr»] = a, (i = 1,2,3) and angular separations, 
#12 = #13 = 27r/3 (see Fig. BJ. As for the trimers, ther- 
mal fluctuations about the ground state configuration 
give rise to significant dipole moments that are absent at 
T = 0. To tackle this issue we estimate the solvation free 
energy of a tetramer by considering the harmonic nor- 
mal modes of angular oscillation about the ground state 
configuration. (Note that these modes already enter into 
the calculation of the corresponding association constant: 
see Appendix El) For each mode, a thermal average of 
the contributions of the individual multipole moments 
is computed; the sum of these mean-square terms then 
provides a value for the overall multipole free energy. Of 
course, this approximation neglects the nonlinear inter- 
actions between the modes but, because of the relatively 
low value of the critical temperature, this should not be 
numerically significant. 

Following our treatment of trimers, we will fix the three 
satellite radii at oi; the exclusion zone for the tetramer 



will be approximated by a sphere of fixed radius 04 which 
we choose to center on the positive core ion (of charge 
+zqo). Ideally, the origin of the effective exclusion sphere 
should again be placed so that, say, the total contribution 
of the quadrupolar free-energy term is minimized. For 
the present calculations, however, only the case z — 3 will 
be utilized: then the tetramers are neutral so that both 
monopole and dipole moments are independent of the 
origin about which they are defined. The variations of the 
quadrupole moments due to small displacements of the 
origin and, likewise, variations in the exclusion diameter 
04 that might reasonably be associated with thermally 
induced shape changes, may be neglected at the level of 
precision appropriate in light of the other approximations 
of the theory. (Note that changes in the definition of 04 
are studied quantitatively in Sec. lVIl below.) 

A planar tetramer has three angular normal modes: 
two 'in-plane' modes and one 'out-of-plane' mode. The 
first two modes correspond to (p — in Fig. Q] [and tp* 
I) in l!A4j) of Appendix E] since the three satellite ions 
remain in the (x, y) plane and ion 1 may be considered 
as fixed on the x-axis (at Xi=r 1 =a 1 ). The first mode 
(a) is a 'flapping' mode in which the ions 2 and 3 (see 
Fig. QJ oscillate in phase, towards and away from the axis 
formed by ion 1 and the central, positive ion; in other 
words one has #12 = Q\ 3 [following from Y = in (Sit ]- 
The second mode (b) is a 'pendulum mode' in which the 
angle between the ions 2 and 3 remains fixed, equal to its 
equilibrium value so that 6*12 + 6>i3 = 47r/3 corresponding 
to X — in JS3J]- Lastly, in the 'out-of-plane' mode 
(c), two satellites are fixed whereas the third one swings 
around the plane of equilibrium (corresponding to the 
mode where X = Y = while tp is varying) pul |. 

For each mode we need a configuration-space weight- 
ing factor: these all derive from the expression 113. 2|) for 
the association constant. For the tetramer, using the co- 
ordinates in Fig. Q] this is 

G?riG?r2<ir3 = 8ir 2 ri 2 dridr2(ir 3 

x {r 2 2 sin l2 d9 12 ){r3 2 sin 9i 3 d9 13 )d<p, (4.28) 

where a prefactor Awri 2 comes from the full angular inte- 
gral over the orientation of the x axis, while a factor 2w 
arises from the axial integral (rotating the y axis about 
the x axis so that satellite 2 is in the (x, y) plane). Insofar 
as we consider the angular modes at fixed ri—ai, the only 
relevant factor for the angular averages over the mode 
coordinates is sin 612 sin 6*13 d&\ 2 d02 3 d(p. (Note that es- 
sentially identical considerations enter in writing H4,25|l 
where 6*12 — 2a.) 

Now the monopole moment of the (general) tetramer 
is always Qo,o = {z — 3)qo/V^', but the higher moments 
clearly depend on the mode configurations as we proceed 
to specify. 

(a) In-plane flapping mode. Let #12 = #13 = 6 be 

the angle describing this normal mode (but recall that 
6 = 2tt/3 specifies the ground state). The dipole and 
quadrupole amplitudes generated by excitation of the 



12 



mode are then 



£ \Qi, m \ 2 = — [1 + 2 cos9] 2 g 2 al, 



4tt 
15 



V \Q 2 , m \ 2 = -r{sm 4 8 + 3cos 4 e} q y i 

' " m. 47T 



(4.29) 
(4.30) 



The reduced repulsive Coulombic energy between the 
three satellite ions is given by 



E a (9) = - — 
zai 



1 



1 



sin(0/2) 2 sin 



(4.31) 



Thus the thermal average square moments may be cal- 
culated from 

(|Q?,ml) = W*) d0sin 2 9 \Q Lm (9)\ 2 e E ^ T ' , 

Jtv/3 

(4.32) 

where l/M a {T*) is the obvious normalizing integral. The 
limits specified on 9 correspond to the hard-core restric- 
tions in the case a\ = a and should be relaxed appro- 
priately if a\ > a (although, since they correspond to 
the maximal interionic repulsions, the differences will be 
small) . 

(b) In-plane pendulum mode: Now let us put 

012 = (2tt/3) - 9' and 9 13 = (2tt/3) + 0', so that 0' de- 
scribes the angular amplitude of the mode. The dipole 
and quadrupole amplitudes are then 



E \Q 

' " m. 



i m| 2 = — [l-COS0']^ a 2 ; 



(4.33) 



3 

2^ 1 

V |Q 2 ,„| 2 = -^[5-2cos(20')] 9o 2 at. (4.34) 

^— <m 167T 



The thermal average is now computed via the normalized 
integration 



(\Qi, m \ 2 )b = M b (T*) f d0'[l + 2cos(20')]x 

J-ir/3 

\Qu m {9')\ 2 e E ^/ T ' , 
where the reduced energy can be written as 



E b {9') = - — 



1 2V3cos(0'/2) 
Vf + 1 + 2 cos 0' 



(4.35) 



(4.36) 



(c) Out-of-plane mode. Finally, in the out-of-plane 
mode as described previously (in which 9i2 = 0i3 = 27r/3 
and <p varies), the dipolar and quadrupolar amplitudes 
are 



J2 \Qhm\ 2 = — {1 ~cos tp}q%al 



8tt 
45 



(4.37) 



V |Q 2 ,m| 2 = — [3-2cos^ + 3cosV]'Z 2 at ) (4-38) 

' — '171 b47T 

while the reduced repulsive energy is simply 



Ec(<p) = - 



1 



C0S(<y9/2) 



(4.39) 



For this mode, the thermal average is performed accord- 
ing to 

(|Q;,™| 2 ) C = N C (T*) f m d V \Qi, m (<p)\ 2 e E °M , (4.40) 



where the condition <p < tp m = it— 2 arcsin(l/v / 3) ex- 
presses the hard-core condition and l/Af c (T*) is the cor- 
responding normalizing integral. 

At this point, the overall solvation free energy of a 
tetramer, f^(T; {p a }), may be calculated by summing 
the solvation free energies computed for each mode. This 
completes the basic general analysis. 



V. CRITIC ALITY AND COEXISTENCE 
UNDER CHARGE ASYMMETRY 

A. Pure Debye-Hiickel Theory 

The original Debye-Hiickel theory Q amounts to writ- 
ing the overall free energy density as 

fir, P+ , P _) = / DH (T, P+ , P _) + r (r, P+ ) + r (r, P _) , 

(5.1) 

where p+ and p_ are the densities of cations (with charge 
<7+ = zqo) and anions (with charge g_ = —qa), respectively, 
and / DH was obtained in l|4.16H • From the electroneutral- 
ity condition 12. lit one has 

p+ = p/(l + z), p-=zp/(l + z), (5.2) 

while using the expression Il4..'>ll for Debye length, with 
x=Ka, the normalized density ll 1.211 becomes 



p* = x 2 T*/4tt. 



(5.3) 



The contribution to the overall chemical potential from 
the DH free energy is then 



P DH = -x/2T*(l + x). 



(5.4) 



Taking C+ = C_ = A? in <|2j| where A X (T) is the de 
Broglie wavelength for free ions, the ideal gas contribu- 
tion is merely 



1 



/i Id = \n(x z T*) + In 



1 



1+2! V 1 + Z 



Z Inf Z 



1 + z V 1 + z 



K4S) ^ 



The overall chemical potential is p — p DH + p Id , while the 
reduced pressure follows from Il2..'ill as 

p i = 4ira 3 p = x 2 T* + ln(l + x) - x + \x 2 /(l + x) , (5.6) 

which is the same as (4.6) in I and quite independent of 
z. 

Since the expression for the pressure does not de- 
pend on z and the overall chemical potential is also z- 
independent except for the constant terms in the ideal 



gas form 115, 5|) . the conditions for criticality and phase 
coexistence are identical to those derived in I for the 1:1 
model. The phase coexistence curves are likewise identi- 
cal: see I Fig. 1(a). In summary, the pure Debye-Hiickel 
theory predicts that the critical parameters are indepen- 
dent of z and given by 

T* = 1/16, p* c = 1/64tt, x c = l, 

Z c = p c /p c fc B T c = 16 1n2-ll, (5.7) 

while the numerical values are presented in Table ITTI 

B. DHBjCIHC Theory 

Extending the DHBjDIHC pairing-plus-solvation ap- 
proach for 1:1 electrolytes to z:l electrolytes, we now in- 
clude dimers, trimers, and all further primary clusters up 
to (z + l)-mers, and add their free energies to the overall 
electrostatic free energy to obtain 

f(T; {P.}) = £ [f d (Pu) + f7(T; { Pa })} , (5.8) 

where v = +, — , 2, 3 . . . for positive ions, negative ions, 
dimers, trimers, . . ., respectively. To determine the de- 
gree of association of the free ions into dimers, trimers, 
. . ., we need the association constants K m ^ z (T*) as com- 
puted in Sec. IIIII Then, under chemical equilibrium the 
cluster densities p m for 2 < m < z + 1, satisfy the laws 
of mass action in the form 

Pm = K m - ljZ p+p™- 1 exp[fj% + (m - 1 V- 1 - fC] > ( 5 - 9 ) 
where the excess chemical potentials are given from I by 

pl\T-{p a })^-dr/d Pu \ Tpai . (5.10) 

We dub this extended treatment DHBjCIHC theory for 
"Debye-Hiickel theory supplemented by Bjerrum associ- 
ation into Clusters that are solvated by the Ionic fluid, 
and Hard Cores." 

In order to obtain a canonical equilibrium state of the 
(z + 2)-component fluid, one needs, in addition to elec- 
troneutrality and the z mass action conditions, one extra 
parameter, such as the overall density p or, more conve- 
niently, the reduced Debye variable x — na. Moreover, 
phase coexistence entails the conditions 112.411 and (22J 5 
whereby one can show that the equality of all the differ- 
ent electrochemical potentials between coexisting phases 
can be replaced by the equality of the chemical potential 
of the neutral species alone (dimer, trimer, or tetramer, 
respectively, for 2 = 1, 2, or 3). (For the charged species, 
the eZecirochemical potential must match between two 
coexisting phases as mentioned in the Introduction and 
discussed in Sec. IVIIll below in connection with the Gal- 
vani potential.) One effective computational strategy is 
thus to plot parametrically (p(x,T), j2 n (x,T)) (where n 




FIG. 7: Examples of the variation of the pressure with the 
chemical potential of the neutral species p„ (shifted by an 
arbitrary constant), calculated for a 3:1 electrolyte treated 
within the DHBjCIHC theory with refined standard param- 
eters. Two-phase coexistence below T c can be realized when 
the curve intersects itself (while, as usual, the states below 
the intersection are not stable). The plots are constructed 
parametrically as functions of x = na using increments of 0.03 
around x c — 1.570 for reduced temperatures T* — 0.04250, 
0.043345 and 0.04380. 



denotes the neutral species) , and to seek for two different 
values of x giving the same point: see Fig. Especial 
care is needed in determining the coexistence curve below 
criticality for z = 3. 

However, calculations in the single-phase region (above 
T c ) are relatively straightforward; consequently, for the 
purpose of calculating T c and p c another useful approach 
is to generate supercritical loci which must intersect at 
the critical point. We choose the maxima of the k- 
susceptibilities (see Sec. IVIIll . which, indeed, lead to fast 
and accurate determinations of T c and p c . 



C. Geometric parameters for the ion clusters 

To proceed further in the quantitative evaluation of the 
electrostatic contributions to the free energy as derived 
analytically in the previous sections, we must address the 
values and thermal variations of the satellite separation 
radii, oi (which should both depend on the cluster species 
to and the valence z) and of the effective exclusion sphere 
diameters a m for to > 2. Both these issues were discussed 
for the basic 1:1 model (the RPM) in I (see Sees. 6.3 and 
7.1) and so will be treated fairly briefly here. 

If we write a\ =a[l + si :?niZ (T)], it is, first, clear that 
si, m ,z vanishes when T — > for all to and z, so that hard- 
core contact is, in fact, rather rapidly approached when 
T falls. Indeed, for T* < 0.055 [~T*(z = l)], the analysis 
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of I indicates that si,i,i decreases almost linearly with 
T* from si,i,i — 0.08. When z>l, because of the tighter 
binding induced by the larger central charges, which is 
only partly offset by repulsions from the remaining m — 1 
satellite ions, one must also have si, m;Z (T*) < Si iTn , )Z /(T*) 
when z' > z. 

Then, one should also observe (see Figs. and Ta- 
ble lUj that for larger values of z (> 2), the critical tem- 
peratures fall, so that the relevant values of si ;TOjZ (T*) 
will again be smaller than for the 1:1 model. Never- 
theless, one must notice from 114. 13j) . 114. 20j) . Q4.24H and 
(I4.29|l . etc., that the dipolar and quadrupolar contribu- 
tions (when the latter do not vanish by choice of the 
parameter p) are proportional to a\ and a\ , respectively. 
However, in compensation, these powers are always ac- 
companied by the inverse powers a~ and a~ 3 , respec- 
tively, of the exclusion diameters a m (m > 2), which are 
proportional to the corresponding values of a\ (T) and so 
act to reduce the overall sensitivity. 

In I, the choice of the radius a,2 for the effective exclu- 
sion sphere that approximates the true bispherical exclu- 
sion zone of a dipolar dimer (see Fig. EJ was discussed 
by considering various bounds and their mean values. It 
was decided to accept, as most appropriate, the 'angu- 
lar average' value, defined as the radius averaged over 
solid angle of the true exclusion zone as measured from 
a symmetrically located origin of a cluster in its ground 
state. For dimers, trimers, and tetramers in their ground 
states, these angular averages are, respectively, 

a% 3 3, n a% 5 a% 11 
-±-=- + -ln3, -2- = -, = — , 

a 4 8 a 4' a 8 

~ 1.16198, = 1.25, = 1.375. (5.11) 



It transpires in the calculations leading to the criti- 
cal parameters, that the predicted values for all three 
cases, z=l, 2, and 3, are dominated by the properties of 
the primary neutral clusters, namely, the neutral dimers, 
trimers, and tetramers, which prove to be by far the most 
abundant species. In turn, for fixed z, these are found 
to be the most sensitive to the geometrical parameters. 
Accordingly, we have examined (as, in fact, did Levin 
and Fisher) various other more-or-less plausible criteria. 
One simple, but clearly rather arbitrary possibility, is to 
choose a a so that the approximating exclusion sphere has 
a volume matching that of the true exclusion zone. We 
identify these parameters as 'steric': they take the values 



i\ _ 3 a\ _ 19 1 / 3 a% _ 7 2 / 3 

~ 1.19055, ~ 1.33420, ~ 1.45220. (5.12) 



Another choice, since the interactions that are being 
truncated by the exclusion zones are Coulombic, is the 
harmonic diameters defined as the inverse of the angular 
average (again taken from the clusters geometric center 
of symmetry) of the inverse radial distance to the surface 
of the exclusion zone. This leads to the values 



a\ _ 6 4 _ 2 a'l _ 4 

~~a ~ (2 + 3 In 3)' ~~a ~ (1 + ln2)' ~~a~ ~ (1 + 3 In 2)' 

~ 1.13297, ~ 1.18123, ~ 1.29894. (5.13) 



Compared to the angular averages i|5.11|l , these are some 
2.5 — 5.6% lower which leads to increased solvation. 
While, in accord with I, we judge that the angular aver- 
ages are to be preferred, the predictions of the steric and 
harmonic parameters will be discussed below. 

In as far as the satellite separation ai(T) exhibits a 
T-dependence, this will be inherited by the a a (T). How- 
ever, in the case of charged asymmetric dimers, as needed 
for z > 2, the offset of the center of the effective exclu- 
sion sphere from the clusters' geometric center [as em- 
bodied in l|4.19fl ] naturally raises the question: Why not 
calculate the a CT 's from the offset center ? Likewise, at fi- 
nite temperature, the crucial bending fluctuations of the 
trimers and tetramers obviously suggest further modifi- 
cations in the calculation of the a CT 's. The temptation to 



explore these refinements however, may be resisted, first, 
because the effects are likely to be small, and, just as 
important, because the resulting changes in critical pa- 
rameter estimates will be less significant than result from 
other approximations already accepted. 



VI. QUANTITATIVE PREDICTIONS 

At this point, it is imperative to re-emphasize that the 
primary aim of the present study is to elucidate the basic 
physical mechanisms underlying the systematic trends in 
the various critical parameters that are induced as z in- 
creases, and, at a semiquantitative level, to understand 
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TABLE II: Predicted critical parameters, T* = Dak B T c / zqg , 
p* = Pcd 3 , x c — K c a, Z c = pc/ pck B T c and the mole fraction 
of free ions, y± c = (N+ + N-)/N\ c , for z:\ hard sphere elec- 
trolytes, as predicted by the DHBjCIHC theory with 'stan- 
dard' parameters (refined for 2 = 3): see text. Monte Carlo 
results fill2Sj are displayed in parentheses. 

z 10 2 T* 10 2 p* c x c Z c y± c 

DH 6.250 0.4974 1 0.9063 1 

1 5.567 (4.93 3 ) 2.614 (7.50) 1.038 0.2451 0.1828 

2 4.907 (4.70) 6.261 (9.3) 1.366 0.1708 0.1164 

3 4.334 (4.10) 11.90 (12.5) 1.570 0.1433 0.0838 



TABLE III: Critical-point mole fractions, y a , of the primary 
clusters (expressed as percentages) according to their total 
charges, q a , for 2:1 models described by DHBjCIHC theory 
(with 'standard' parameter values). Unlabeled clusters are 



monomers. 
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-1 
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+2 
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81.72 


9.14 
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2:1 


10.33 


72.93 


15.43 
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(trimer) 


(dimer) 






3:1 


8.04 


77.17 


11.13 


3.32 


0.34 






(tetramer) 


(trimer) 


(dimer) 





the magnitudes of the changes. Recall that the true val- 
ues of T*(z), etc., are already known to satisfactory ac- 
curacy from the recent simulations 0,E3]- Consequently, 
a uniform theoretical treatment of the 1:1, 2:1, and 3:1 
models is of greater importance than are concerns for 
various specific subtleties that we know, a priori, cannot 
yield truly reliable and accurate critical-point data owing 
to our failure (not to say inability) to treat adequately 
the essential critical fluctuations: see, e.g., [20]. The fluc- 
tuations, of course, serve to realize the universality class 
of the critical behavior 0, Q H 0, 0, El while, at the 
same time, depressing the critical temperature and (for 
these primitive electrolyte models) increasing the criti- 
cal density relative to the predictions of even "the best," 
classical mean-field, or self-consistent treatments. 

With these points in mind, the principal explicit nu- 
merical calculations of the electrostatic free energy terms 
that we have undertaken have utilized the simple (T = 0) 
angular averages diameters a a , listed in H5.11I) and, fur- 
thermore, have accepted the "in-contact" or T — ► limit, 
a\ = a, for the satellite ion separations in all clusters. 
It should be stressed, however, that the calculations of 
the cluster association constants, K m<z (T) in Sec. lIIII are 
not so constrained: rather, each satellite ion is allowed 
to explore the full phase space restricted only, at large 
separations, by the Bjerrum-type optimal cutoffs, i? m , z . 

A. 1:1 or Restricted Primitive Model Electrolyte 

Here, we merely refine the results of Fisher and Levin 
from I. Within the DHBjCI theory using the angular av- 
erage for <X2 (but without the contribution / HC ), one finds 

T* = 0.05740 and p* c = 0.02779. (6.1) 

We may supplement the results of I by recording that the 
use of the larger steric parameter af [see H5.12I) ] modi- 
fies the predictions for T* and p* by factors 0.9815 and 
1.0086, respectively, whereas, the smaller harmonic ex- 
clusion diameter aty, yields factors 1.0195 and 0.9989. 
Hence, a larger cluster size leads naturally to a decrease 
in T*, since fewer attractions are realized, and to an in- 
crease in p c . 



The next step is to include the hard-core term / HC . 
Keeping the angular average radius a?; and taking the 
bec hard-core value £? -/a^ = 4/3"v / 3, a choice of parame- 
ters that we will refer to as 'standard', we find the critical 
parameters displayed in Table ITTI As expected, the intro- 
duction of hard-cores reduces both the critical temper- 
ature (by around 3%) and the critical density (by 6%). 
These effects are stronger if the low-density limiting value 
B CT /a^ = 27r/3 is used since T* then drops to 0.05293, i.e. 
by 5%, whereas p* becomes 0.02469, falling by 6%. Fi- 
nally, using the angular average a§ but with the choice 
B a ja? a = 1.300, which lies between the low-density and 
bec values, we obtain the 'optimal-fit' estimates 

T*(z = l) = 0.05455 and p* c {z = 1) = 0.02542 . (6.2) 

The coexistence curves predicted by the DHBjCI the- 
ory (with the angular average for a 2 ) and by the DHB- 
jCIHC theory with standard parameters are plotted in 
Fig. 03 The introduction of / HC significantly lowers the 
liquid sides of the coexistence curves. One may notice 
that the Monte Carlo data could well be better fitted by 
some choice of B a between and the bec value. 

Note that in addition to the 'standard' values of T* 
and p*, listed in Table ITTI the last column, labeled y± c , 
reports the critical value of the mole fraction of unasso- 
ciated ions, namely, 

y±=JJ+ +V-, with y <y ^n a N a /N , (6.3) 

where y a is the mole fraction of species a while n a is 
its ionic weight (i.e., n a = 1 for free ions but n a = m 
for a cluster of one positive charge and (m — 1) negative 
charges). For the 1:1 model, we have y+ c = y- c = 0.0914 
while the critical mole fraction of the associated ion pairs 
is y 2c = 2/?2 c/p = 0.8172: see Table ITTT1 The fact that 
(within the DHBjCIHC theory) almost 80% of the ions 
are associated into dipolar ion pairs near criticality makes 
it less surprising that a model of neutral but charged hard 
dumb-bells might have a comparable coexistence curve, 
as some simulations suggest [45, 46]. 

Within this DHBjCIHC approach, we remark that an 
increase in the value of the association constant K11 
yields a decrease of T* for a 1:1 electrolyte: with the 
standard parameters, we find that varying K\ t i by ±5% 
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FIG. 8: Coexistence curves computed for 1:1, 2:1, and 3:1 eq- 
uisize hard sphere ionic fluids or primitive model electrolytes: 
the solid lines correspond to the DHBjCI theory (without ex- 
plicit hard-core excluded- volume terms); the dashed curves 
include 'standard' bcc hard-core terms. The exclusion di- 
ameters used are the angular averages l|5.11fl except for the 
refinement 04/a = 1.41 for the 3:1 model. Solid symbols rep- 
resent simulation estimates for the critical points |7l \2h\ and 
open symbols the coexistence curves based on precise RPM 
simulations Q. 



around the value l-i.ot results in changes of T*(z = 1) 
of order =F-0.2%. While changing an association constant 
does not affect directly the free energy of our model (in- 
deed, T* is the same in DH theories with or without ion 
association [13), it nevertheless affects the various mole 
fractions, which do enter in the solvation free energies 
f™. To our knowledge, the variation of the sum J^a fa 1 
on increasing the association constant can only be deter- 
mined post-facto: for z — 1, our results indicate a more 
weakly coupled system with a lower critical temperature, 
in accord with the findings of Jiang et al [ij; however, 
for z = 2 and 3, we find the opposite trend on varying 
K z , z , as noted below. 

Let us also recall that I presented numerical and graph- 
ical data showing how the predicted values of the critical 
parameters depend on the choices made for the mean 
ion separation a\ and for the exclusion radius These 
results may reasonably be taken as indicative of the cor- 
responding shifts that are likely to arise in our analysis 
of the 2:1 and 3:1 models. 



B. 2:1 Hard Sphere Electrolyte 

We report first the basic DHBjCI results, using the 
angular averages a 2 and 03 listed in H5.11J1 : they are 

T* (z=2) = 0.05235, and p* c (z = 2) = 0.06429 . (6.4) 



The corresponding coexistence curve is plotted in Fig. |H1 
One might note, first, that as in the 1:1 model, the shape 
of the liquid side of the coexistence curve below about 
0.9T c (z = 2) becomes markedly concave. This behavior, 
while violating no known thermodynamic or other con- 
ditions, certainly appears unphysical. Furthermore, by 
comparison with the true results indicated by the simu- 
lations, this concavity must be judged as quite mislead- 
ing. No doubt it results from the failure to satisfactorily 
describe the correlations, and thence, the free energy of 
the low-temperature liquid at densities p* > 0.15 via a 
collection of free ions plus fairly compact neutral and 
singly charged clusters. This general issue is also ad- 
dressed briefly in I Sec. 8.5: it may be noted that the 
standard MSA exhibits similar although somewhat less 
pronounced features: see I Fig. 8(d). 

On the other hand, the downward shift in T c and the 
marked increase in p c reproduce most satisfactorily both 
the trends and the magnitudes obtained in the simula- 
tions 28): see Figs. ^ and El These trends are also 
reproduced fully by the other choices of exclusion diam- 
eters. However, as could be expected, the sensitivity to 
the size of the bigger clusters is enhanced in the 2:1 case 
compared to the 1:1 model. Indeed, on using the steric 
diameters, we find T* — 0.04850 and p* — 0.0737 imply- 
ing a drop by 7.3% and an increase by 15%, respectively. 
With the harmonic parameters, the conclusions are re- 
versed yielding T* = 0.05574 and p* c = 0.06012. As dis- 
cussed below (and see Table |Hj , a larger fraction of the 
ions are bound in the clusters when z = 2 compared to 
z = 1, thereby amplifying the sensitivity to the cluster 
characteristics. 

This enhanced sensitivity is also found for the hard- 
core effects: thus the values of T* and p% predicted by 
the standard DHBjCIHC theory, listed in Table HT| are 
lower by 6% and 3% relative to the values in <6.4ll . Like- 
wise, the low density value of B a yields T* =0.04375 and 
p* = 0.06422, which seriously overestimates the hard-core 
effects, making T* drop to well below the Monte-Carlo es- 
timate. However, the 'optimal-fit' choice B a /a^ = 1.300 
yields the values 

T* (z = 2) = 0.04691 and p* c (z = 2) = 0.06285 , (6.5) 

which, indeed, provide the best fit of our analysis to the 
Monte Carlo data (see Figs. and EJ. However, the cor- 
responding coexistence curve is excessively narrow even 
compared to the standard prediction shown (dashed) in 
Fig. 

On the other hand, it transpires that the sensitivity 
to the association constant is not so great. Thus in 
the DHBjCI approach, changing the cut-off for i^2,2 by 
±20%, induces changes in ^2,2 of order ±1%, leading to 
shifts in T* of order ±0.003%' and in p* of order ±0.2%, 
totally negligible within our level of approximation. 

Returning to the standard DHBjCIHC theory, one sees 
from Table [H] that it predicts a drop in T* (compared 
to the 1:1 electrolyte) of order 12% and an increase in 
p* of 140%. These results are to be compared with the 
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Monte-Carlo results indicating a drop in T* of 5% and an 
increase in p* of around 24%. The predicted T* and p* c 
agree within 4% and 33%, respectively, with the current 
Monte Carlo estimates. The overall quantitative results 
are therefore fairly close to the Monte Carlo values, indi- 
cating that the main physical features have been captured 
by the theory. 

As regards the composition of the fluid at criticality, 
one learns from the last column of Table [H] that fewer 
than 12% of the ions now remain free or unassociated, 
even less than predicted in the 1:1 case. As can be seen 
from Table HHl the free +2go ions are strongly depleted, 
less numerous than the —qo anions, by a factor 1/8. In- 
deed, the numbers of positively charged dimers roughly 
match the oppositely charged free anions. However, while 
the predicted overall association rate is larger than for 
the RPM, the fraction of the ions bound into the neu- 
tral or "molecular" clusters (now trimers) is some 11% 
smaller. Needless to say, the values of the y a listed in 
Table UTA verify the electroneutrality condition that im- 
plies 2y + + ±y 2 -y- = 0. 

C. 3:1 Hard Sphere Electrolyte 

As before, let us first record the predictions of the basic 
DHBjCI theory using the angular averages, needed now 
for 02, 03, and 04, the last for the tetramer which we 
expect to be the dominant species near criticality. We 
find 

T* (2 = 3) = 0.05054 and p* c (z = 3) = 0.1063 , (6.6) 

where the corresponding coexistence curve is again dis- 
played in Fig.|Hl These values, as is also clear from Figs.Q] 
and[3 continue to reproduce the appropriate trends with 
increasing z as originally revealed by the Monte Carlo 
simulations. However, as also evident in Fig.^ the drop 
in T* of only 3.5%, relative to the 2:1 model, is signifi- 
cantly less than indicated by the simulations: in fact, the 
result l|6.6Jl suggests a concave variation for T*(z) rather 
than the convex behavior maintained by the simulations 
up to z = 4 HHH]. 

The 'culprit' is obviously the failure to take explicit ac- 
count of the hard-core excluded volume effects important 
above and even at the predicted critical density which is 
65% larger than for the 2:1 model. The very slow decay 
of the liquid side of the coexistence curve when p in- 
creases, as seen in Fig. [SI strengthens the point. Indeed, 
the standard DHBjCIHC theory yields 

T* (2 = 3) = 0.04580 and p*(z = 3) = 0.1089 . (6.7) 

Relative to the 2:1 model the critical temperature has 
now fallen by 12.5%, which may be compared with the 
Monte Carlo drop of 12.8% (see Table However, the 
value of p c has changed rather little. 

These predictions are not quite those entered in Ta- 
ble |n] because it was deemed worthwhile for this case to 



explore further the influence of the exclusion radii. In 
particular, as discussed in I Sec. 6.3, the mean size of a 
physical cluster, for any sensible definition will grow with 
increasing temperature. Hence, in choosing the tetramer 
exclusion radius 04, it is reasonable to consider for use 
near T c a value somewhat larger than the T — angular 
average of = 1.375a [see H5.11I) ]. Having examined the 
effects on the values of both T* and p* , the ratio 

o 4 /o = 1.41Q, (6.8) 

was selected as a preferred refinement of the standard 
parameters. (The increase of 2.6% brings the ratio to 
almost midway between a\ja and the steric value a|/o~ 
1.452.) Accordingly, Q6.8J1 has been adopted for comput- 
ing the results displayed in Table ITT1 in Figs. and and 
elsewhere below; the corresponding coexistence curve for 
z = 3 is displayed (dashed) in Fig. El 

Evidently, the trends observed as z increased from 1 to 
2 are now continued regularly; and the previous concave 
variation of T c (z) is no longer so apparent. Furthermore, 
the trends still mirror rather faithfully those given by the 
simulations: These indicate a rise in p* by 35% when z 
changes from 2 to 3; the standard calculations yield a 90% 
relative rise which is significantly greater, but, as seen in 
Fig. 121 not at all unreasonable. Indeed, T* agrees with 
the simulations to within 6% while p* agrees to within 
5%. Overall, both the magnitudes of T*(z) and p*(z) 
and the trends with z must be judged quite successful! 

From Table IIIII we see that the overall fraction of free 
ions remaining at criticality has now dropped still further 
to about 8.4%. At the same time close to three quarters 
of the ions are again bound in the neutral, molecular clus- 
ters (now tetramers). The fraction of free, unassociated 
cations of charge +zq n continues to fall dramatically as 
z increases: on a heuristic basis, a decay like y+, c ^e~ bz 
seems not implausible. Following the thought of Shelley 
and Patey |45j, one might also speculate that a system 
of rigid, neutral molecules or (z+l)-mers formed of z+1 
equisize hard spheres with z of charge —qo attached sym- 
metrically to a central sphere of charge +zqo, might con- 
tinue to mimic the z:l equisize hard-sphere ionic systems, 
at least up to z < 12. Beyond that, packing effects in the 
satellite ions could play an important role. 

It is probably appropriate to point out, as anticipated, 
that our z = 3 predictions are less robust than those for 
z < 2. Thus the 'optimal-fit' assig nment B a /al = 1.300 
together with the angular averages a?; and Og but taking 
a4/a = 1.390 yields 

T* (z = 3) = 0.04136 and p* (z = 3) = 0.1171 , (6.9) 

which reproduces the Monte Carlo results quite satisfac- 
torily. However, this choice once more leads to a coexis- 
tence curve which falls much too steeply when T < 0.9 T c . 
Again, the low-density value for B a gives T* < 0.038 well 
below the simulation value. 

On the other hand, if one uses the [5/2] Pade approx- 
imant for the association constant integral in ij3,12|) 
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in place of the more accurate fit ij3,13ll one finds that the 
resulting 4% decrease in K 3j3 [see Fig. E] leads to a de- 
crease in T* of only 0.14%. The effect on p* is similar 
and hence, post facto, of little consequence. 

Finally, it is interesting to note from Table ITT1 that the 
Debye length at criticality, namely, £ DiC = 1/ ' K c = a/x C) de- 
creases steadily as z rises. In essence, this merely tells us 
that larger central charges in ionic clusters lead to tighter 
screening. It should be noted, however, that £ D (T, p), as 
defined in 114.311 is not really susceptible to either physical 
measurement or simulation since our definition depends 
on having a well defined, but intrinsically somewhat arbi- 
trary decomposition of the system into distinct species of 
ionic clusters. On the other hand, the prediction that the 
critical pressure ratio, Z c =p c / p c k B T c , decreases strongly 
as z increases (see Table ITTI column 5) should be open to 
test by simulations. 



VII. SPECIAL INFLECTION LOCI 

In determining numerical values of critical parameters 
from a given model free energy, it is natural to start by 
calculating the two sides of the coexistence curve, pi(T) 
and p v {T), using techniques, such as illustrated in Fig. [3 
in principle, one can then raise T and monitor Ap = 
pi (T) — p v (T) , determining T c from the vanishing of, say, 
Ap 2 and, then, p c from, say, \(pi — p v ) evaluated at 
T c ost . In practice, however, this method proves tedious 
and as experience (and a consideration of Fig. EJ reveals 
is poorly adapted for providing precise and accurate (i.e., 
reliable!) values of T c and p c . 

An effective alternative is to confine attention to the 
one-phase region above T c where, in the first place, calcu- 
lations are more straightforward since, in particular, no 
'two-phase solutions' need be sought. Then, as demon- 
strated recently in simulations [fll.l3ll.l32T] (although also 
of value in studying experimental data) one may seek 
various loci, say p^(T), which all converge on the crit- 
ical point. Since the isothermal compressibility \t — 
(d p I 8p)t I p diverges at criticality, one obvious such locus 
is provided by those densities, say po(T), on which, at a 
fixed temperature above T c , the compressibility achieves 
its maximum. But by considering the inflection points of 
the standard isothermal plots of p vs volume or vs den- 
sity, one soon realizes that this locus is but one of a nat- 
ural family of fc-loci, say p (k) (T), HJUllIzL on which 
the so-called fc-susceptibilities, x P) = x(T,p)/p k , 
attain their maxima: equivalently, these are just the loci 
of isothermal inflection points of plots of p vs p k . 

Because of their potential usefulness in simulation and 
experiment, the behavior of the fc-loci in the scaling re- 
gion close to criticality has been investigated in some 
detail I n the case of general, nonclassical criti- 

cal points, they exhibit nontrivial and informative singu- 
lar behavior as functions of t = (T — T c )/T c as fc varies. 
However, for classical critical behavior, as relevant here, 
all the fc-loci asymptotically approach the critical point 



(T c , p c ) linearly in the (T, p) plane. Thus by numerically 
determining two or three loci — for the results reported 
here we utilized fc = 1, 0, and — 1 — and solving for their 
mutual intersection point, one may locate T c and p c . In 
practice the method proves efficient and precise. 

More generally, however, the nature of the loci further 
from criticality and a possibly characteristic dependence 
on z is a matter of interest to which we now turn. 



A. Debye-Hiickel predictions 

To gain a little perspective, let us examine, first, pure 
DH theory (as presented in Sec. IV A|l where analytical 
calculations are feasible. Three cases arise as illustrated 
in Fig- EI When fc=l, the pressure isotherm always has 
an inflection point above T c , characterized by na = x — l. 
The (fc = 1) locus is thus a straight line starting at 
(T c ,p c ), namely, p^ *(T) = T*/47r. When fc > 1 one sees 
that by construction, \^ diverges to +00 when p — > at 
fixed T; but this divergence competes with the localized 
maximum driven by criticality. As a consequence, for T 
above but not too far from T c , when p drops beneath p c 
one first encounters a maximum in \ an d then a min- 
imum before the divergence as p — > 0. However, as T is 
raised at fixed fc one eventually encounters an annihila- 
tion or terminal point (T a ,k,p a ,k) at which the minimum 
and maximum merge and the fc-locus is terminated with 
a horizontal slope, i.e. a tangent parallel to the p axis. 
Above T a fe the susceptibility x^ fe ^(T, p) falls monotoni- 
cally as p increases and no 'critical maxima' are realized. 

When fc < 1, a similar scenario emerges for p > p c . 
As a result there is, overall, a termination boundary in 
the (p, T) plane with a minimum at the critical point. 
For DH theory this takes the form of the bold curve in 
Fig. El For fc < 1 the termination boundary approaches 
asymptotically the line T* k(<1 s ) « Airp*/(2 + \/3) 2 while 

for fc > 1 and large p one has T* fe ( >X ) ~ Anp* / (2 — a/3) 2 . 
As also clear from Fig. for values of fc differing much 
from 1, the fc-loci are rather short (and hard to locate 
numerically). The asymptotic slope of the general fc-locus 
at criticality is given by 

(dpW/dT) c = (2/7r)(fc -fc), (7.1) 

where, within DH theory one has fco = 9/8 (for all z). 

B. DHBjCIHC Predictions 

How are these fc-loci affected when the DH approxima- 
tion is supplemented by association, solvation and hard- 
core effects, and how do they evolve with zl Some results 
obtained with the full theory are displayed in Fig.EH As 
expected from our analysis of the DH theory, most of 
the fc-loci do indeed terminate with a horizontal slope at 
some point within the range of investigation. Most of the 
values of fc examined are smaller than 1 and the fc-loci 
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FIG. 11: Plots of the reduced interphase Galvani potential 
A(f> = qoA(j)/kBT for az:l electrolyte as predicted by the pure 
Debye-Hiickel theory. 
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(8.3) 



By substituting in Ij8.ll) and solving for the electrostatic 
potential difference, we obtain 



z In 



P-i 

P-v 



■In 



(8.4) 



On us- 



where A(j)=q Q Acj)(T)/k B T and A(f> = <p Uq - <p 
ing the electroneutrality constraint, we obtain the much 
simpler form 



A$(T) = (l-z- 1 )ln[p l (T)/p v (T)] 



(8.5) 



As anticipated, the predicted Galvani potential A<fi van- 
ishes identically when z = l. 

Fig. El presents plots of this Debye-Hiickel result for 
A(f> vs. T* for various values of z. Note that when T 
approaches the form 118. 5J1 implies that A<j>(T) should 
approach a constant value since p v (T) vanishes exponen- 
tially fast with 1/T Q. We should remark that within 
DH theory the ratio pi/p v is independent of z at fixed 
T. By expanding pi(T) and p v (T) around p c in powers 
of t = (T c — T)/T c , one finds A(j) m , with since the 
theory is classical, j3=h. 

In this simple DH analysis, the Galvani potential is 
rather trivially proportional to the logarithm of the ra- 
tio of densities in the two coexisting phases. One might, 
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FIG. 12: The reduced interphase Galvani potential plotted 
vs T/T c - The solid lines show the predictions of the DHBjCI 
theory (with B a = 0), the dashed curves, the DHBjCIHC 
theory (with refined standard parameters), and the dotted 
plots, the DH theory. 



perhaps, suspect that this indicates the existence of some 
simple "universal" result not depending significantly on 
the detailed microscopic interactions. For a better un- 
derstanding sho wing that this idea is false, let us, fol- 
lowing Bjerrum allow for the formation of dimers 
by association, neglect all solvation effects arising from 
their dipole and higher moments, and treat the dimers 
as of the same size as the free ions. This allows us to use 
the standard DH free energy H4.16J1 for the electrostatic 
contributions, and thence to write the total free energy 
density as 



/(T*;{p CT }) = / DH (T*;{ Pff })+ £ /«( Po 



(8.6) 



Now the electrochemical equilibrium conditions always 
apply and thus, Eqs. I|8.1|) give the Galvani potential cor- 
rectly. Notice, however, that for z > 1, the dimeric ion 
pairs carry a net charge (z — l)qo so that, although elec- 
troneutrality must still be respected in both phases, the 
simple ratio p+/p-, will, in general, be different in the 
liquid and the vapor. Consequently, the simple result 
118. 5JI no longer applies! Clearly, the ratio of p+ to p_ 
depends on the density, pi (T) , of the dimers in the two 
phases. This, in turn, must depend via the mass-action 
laws, on the association constant Ki >z (T) of the dimers 
which then determines the overall degree of association, 
say, ct2 7 (T) , which will vary very differently in each phase 
7. Accordingly, we may write P2 = ct2{T)p and impose 
electroneutrality in both phases to simplify 118. 4J1 . The 
result for A(j> may be written 
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which, by comparison, demonstrates that, in general, 
the simple form 118. 5JI must be modified by nontrivial 
temperature-dependent terms that depend on the de- 
tails of the ionic interaction, etc. Nevertheless, the pre- 
dicted leading temperature variation will still reflect the 
t@ form characterizing the coexistence curve. The anal- 
ysis leading to Ij8.7l) involved only the formation of non- 
neutral dimers; but it is clear that in any realistic treat- 
ment there will be a variety of charged species present 
in temperature-varying proportions determined by mi- 
croscopic details. Thus simple results for the interphase 
Galvani potential should not be anticipated. 

On the other hand, from our explicit calculations of the 
coexistence curves for the 2:1 and 3:1 models, we may de- 
termine A(j>(T) via Ij8.ll) . merely by computing the differ- 
ence in /i_ in the vapor and liquid phases (which quantity 
arises naturally in the computations). The results are 
presented in Fig. El where the temperatures have been 
normalized by the respective critical temperatures to fa- 
cilitate comparison. The plots are qualitatively similar 
to those predicted by the pure DH theory. However, we 
note that in the full theory with standard parameters, it 
is not possible to draw a conclusion regarding the trend 
of Atfi with z. 

The observability of A<j>(T) in a real system is elusive 
if not in principle impossible however, it seems 

that it should be possible to measure A</>(T) in simula- 
tions. Specifically, the potential distribution theorem of 
Widom |4<| OhJ provides a direct way of sampling the (ab- 
solute) electrochemical potential via a suitably weighted 
average interaction of a "ghost test particle" with the in- 
teracting ions in the system which do not "see" the ghost 
particle. The electrochemical potential of a (ghost) ion 
of specific charge should thereby be open to estimation 
in liquid-like and vapor-like simulations of the restricted 
primitive models (or more general models). The appro- 
priate difference should then provide a value of A0(T) . 



IX. DISCUSSION 

Our aim has been to understand, both qualitatively 
and semi-quantitatively, the role of charge asymmetry 
in the criticality of electrolytes. We have extended the 
DHBjDIHC theory of Fisher and Levin H G3 for 1:1 
electrolytes to 2:1 and 3:1 electrolytes by accounting for 
association of ions into charged clusters and including 
the interaction of the clusters with the screening ions 
(solvation). Thus we have labeled the extended theory: 
DHBjCIHC, where the CI now stands for the cluster- 
ion interactions and, for a z:l system, explicit account 



has been taken of the monomers, with charges — qo and 
+zqo, of dimers, trimers, . . ., up to neutral (z + l)-mers. 
The principal results, summarized in Figs. and EJ in- 
dicate that the reduced critical temperature, T*{z), de- 
creases while the critical density increases with increas- 
ing charge asymmetry. Furthermore, these trends and 
the magnitudes of the changes with z agree with the be- 
havior revealed by computer simulations and present a 
significant improvement over the original DH theory. 

To understand the results in physical terms, consider, 
first, the pure DH theory which predicts that the critical 
temperature and density are independent of charge asym- 
metry: as shown in Figs.njandEJ The only direct attrac- 
tive interactions accounted for in this theory are those 
between the ions of opposite charge. These induce a De- 
bye screening cloud around each (monomeric) ion and the 
associated 'solvation free energy' drives the vapor-liquid 
phase separation below T c (z), the vapor phase being sta- 
bilized by the greater entropy available at low densities. 

The temperature is appropriately normalized by the 
energy of attraction of the opposite ions at contact, 
namely, e= \q+q-\/ Da. Under this normalization [see 
(ll.lfl ] the maximum strength of the attractive interac- 
tions is always e and the pure DH theory therefore pre- 
dicts that the (reduced) critical temperature, T*(z), is 
independent of z. 

The DHBjCIHC theory, however, also takes into ac- 
count the formation of ion clusters and treats them as 
distinct species, albeit in mutual chemical equilibrium 
which calls for the calculation of association constants. 
For example, for 2:1 electrolytes, the dimers are species 
with charge +qo while the trimers are neutral. The two 
principal attractive interactions in this case are of mag- 
nitude e between the positive and negative free ions, but 
only ~ \e between the dimer and the negative ion. (The 
interaction magnitude is not precisely \e since the dimer 
has a different effective exclusion zone radius, i.e., a 2 7^ a.) 
Thus, relative to the 1:1 case, the effective attractions are 
smaller for 2:1 electrolytes which explains why the critical 
temperature should be expected to decrease. The same 
argument applies for larger z when there are more inter- 
mediate positively charged species between the free posi- 
tive ions and the neutral clusters. The strongest interac- 
tions is between two free, oppositely charged monomers 
and is always of magnitude e: thus the overall effective 
interaction decreases with increasing z and the critical 
temperature decreases correspondingly. 

To understand the trend exhibited by the critical den- 
sity, p* c {z) = p c a 3 , one must focus on the role played 
by the neutral clusters. As originally shown by Fisher 
and Levin for 1:1 electrolytes, the association of free ions 



into neutral dimers is highly significant at criticality. In- 
deed, according to our theoretical estimates they consti- 
tute about 82% of the overall ion density. For 2:1 elec- 
trolytes, our analysis likewise indicates that about 73% 
are bound in neutral trimers while for 3:1 systems the fig- 
ure is 77% for the neutral tetramers. As a consequence, 
not only are the relative effective charges of the charged 
species decreased by association (as just argued) but, in 
addition, the overall effective fraction of ions in charged 
clusters is diminished when z increases. In leading ap- 
proximation the solvation of a given cluster (charged or 
neutral) is achieved only by charged species. To obtain 
comparable solvation free energy therefore necessitates 
higher overall ion densities and, thereby, an associated 
increase in critical density. The effect is reflected more 
concretely in the expression 114.31 for the effective inverse 
Debye length, k(T, {p a }), whose critical value similarly 
increases with z: see Ta.ble ITT1 

This accounts for the trends displayed in Fig. |2 It 
is interesting to notice, however, that while the Monte 
Carlo results display similar increases in p*(z), the mag- 
nitudes of the increases are rather smaller. It seems likely 
that this is associated with our neglect of the solvating 
influence of the neutral clusters which may be envisaged 
as contributing to a change in effective dielectric con- 
stant. However, the increasing sensitivity of the results 
to the explicit hard-core excluded-volume terms when z 
increases must also be noted. 

It is also appropriate to recall here that our present 
analysis takes no account of critical fluctuations. Exten- 
sive studies demonstrate that the effect of the fluctua- 
tions is to lower T c by 5-10% or more relative to basic 
mean-field-type theories while having little effect on the 
slope of the coexistence curve diameter. In addition, the 
coexistence curve is flattened (since < k). As evident 
from Fig. H the present calculations are quite consistent 
with these general expectations. 

Our results for T c (z) and p c {z) are contrasted with 
those of other available theories in Figs.E3 an dEl For 
this comparison, we have used the standard parameters 
of the DHBjCIHC theory (with a refinement for z = 3) 
as described in Table Inland Sec. lVIl We may note, first, 
that the mean spherical approximation (MSA) |18l llfij. 
like the original DH theory [§], predicts that T* and 
p* remain independent of z. This seems primarily due 
to the failure to take ion association in a sufficiently 
explicit way. A field-theoretic expansion approach ad- 
vanced by Netz and Orland (NO, dashed curves) [2<j], in 
which the particle hard-cores are represented by a sharp, 
large- wavelength cut-off, predicts that T*(z) increases 
strongly with z while p* a (z) falls precipitously at small 
z (< 1) and then rises slowly. In fact, the only previous 
theory known to us that matches the sign of the trends 
revealed by the simulations is the symmetric Poisson- 
Boltzmann (SPB, crosses) integral equation analyses by 
Sabir, Bhuiyan, and Outhwaite jl8|. However, not only 
is the critical temperature for the RPM predicted by the 
SPB theory significantly too high (at T* = 0.0715) but the 
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FIG. 13: Reduced critical temperature, T*(z), as a func- 
tion of the charge asymmetry parameter w = (z — l)/(z + 1) 
as found by simulations (open circles) (2 71, I2«| compared with 
the present calculations (filled circles and triangles) and other 
current theories: MSA SPB, and MPB 00], a 'new mean- 
field theory' (NMF) [33, and a field theoretic expansion (NO, 
dashed curve) [23; note that the predictions of the field theo- 
retic approach have been divided by 10 to bring them within 
the compass of the figure. 



proportionate changes with z are quantitatively much too 
small (by factors of 5.6 and 6.4 for the 2:1 and 3:1 mod- 
els, respectively). Furthermore, the modified Poisson- 
Boltzmann (MPB) approach developed by the same au- 
thors, which they argue should be quantitatively and 
qualitatively better than the SPB, predicts the opposite 
trend for T*(z). Finally, we note that recently devised 
mean field theories based on Kac-Siegert-Stratonovich- 
Hubbard-Edwards transformations of the Boltzmann fac- 
tor [30j , lead to critical temperatures which increase sig- 
nificantly with charge asymmetry again in strong con- 
tradiction to the Monte Carlo estimates (open circles in 
Figs. [H and EH. 

While our theoretical analyses have been based upon 
fundamental principles and provide insight into the vari- 
ation of the critical parameters of charge-asymmetric 
primitive model electrolytes, it must be recognized that 
the results rest upon various approximations. Thus, one 
of our main approximations entails the choice of an equiv- 
alent sphere to represent the exclusion domain of a clus- 
ter. Moreover, we have not explicitly considered higher- 
order association. That, despite these and other approx- 
imations, we find both the correct trends and reasonable 
quantitative agreements with the Monte Carlo simula- 
tions, reinforces our conclusion that the main physical 
features linked to charge asymmetry have been appropri- 
ately captured by the theory. 
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FIG. 14: Simulation estimates (open circles) for the critical 
density, pt{z), compared with those of the present calculations 
(filled circles) and of other approaches: the labels, symbols, 
etc., have the same significance as in Fig. 1131 
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APPENDIX A: ASSOCIATION CONSTANT FOR 
THE TETRAMER 

Consider a cation with charge q+ — zqo at the origin. 
Without loss of generality let the first satellite charge — q$ 
be on the x-axis at ri = (n, 0, 0) in Cartesian coordinates 
as shown in Fig. 0] Taking advantage of the azimuthal 
symmetry, let the second satellite charge be in the x-y 
plane at r 2 = (ra cos Q\2 , r 2 sin 612 , 0) , where 812 is the an- 
gle subtended by the satellite pair (1,2) at the origin. 
Then the most general coordinates for the third satel- 
lite are r 3 = (r 3 cos #i 3 , — r 3 sin #i 3 cos ip, — r 3 sin #i 3 sin tp) , 
where #i 3 is the angle subtended by the satellite pair 
(1,3) at the origin and <p is the angle between the (1, 2) 
and (1,3) planes with ip = representing the planar con- 



figuration. 

The ground state is clearly given by n = r 2 = r 3 = a, 
612 = &13 = 27r/3, and tp — Q. Noting that the main contri- 
bution to the integral defining -ftT 3;3 [see 13.21 1 comes from 
near T — 0, it is helpful to define rescaled coordinates 

0* =(0 13 -27r/3)/\/2T% (Al) 
k={ n /a-l)/T\ 

for i — 1, 2, 3. Then, by expanding about the ground state 
configuration for small T*, one can write the configura- 
tional energy to leading order as 

E3.Z _ 3C 3;Z ^ ^ 7 _ \ * 2 
T* ~ T* ° 3 < 2 Z^ gyg^ 

-^^(ef + ef + eie* 2 ) + o(Vf^), (A2) 

where C 3j2 = 1 — I/Viz. The infinitesimal phase-space 
volume can likewise be written 

dridr 2 dr 3 = a 9 T* 3 dhdl 2 dl 3 

x 8Tr 2 sm 2 (2ir/3)(zT*) 3 / 2 dd*dd;d<p*[l + 0(T*)} . (A3) 

To evaluate the defining integral in ll.3.2li we diagonalize 
the angular quadratic form in (|A2j) by introducing coor- 
dinates 

X={dl+0l)/V2 and F= (6»2 - 0D/V2 , (A4) 
to obtain 

ef + ef + 010; = \{ c ix 2 + Y 2 ). (A5) 

The integrals in l|3.2H can then be evaluated in the form 
H3.3JI . with Jacobean and eigenvalues 

j3 ^< and {A3 < fe M^'2i^^}- (A6) 

Interestingly, the ^T* corrections (and all subsequent 
half-integer power-law corrections) arising in l|A2J) vanish 
upon integration because they are all associated with odd 
powers of the angular variables. An expansion for X 3j2 
can then be found by carrying the expansions in l!A2j) and 
(IA3|I to higher order in T* , and, hence, to higher orders 
in the k and in ip, 9\ and 62- The resulting Gaussian 
integrals can be performed — analytically in low orders 
and numerically, with increasing difficulty, in the higher 
orders — leading to the asymptotic expansion H3.12I) . 

APPENDIX B: MONTE CARLO EVALUATION 
OF THE TETRAMER ASSOCIATION 
CONSTANT 

To evaluate K 3 3 numerically, in order to validate the 
Pade approximants constructed from 113. 1211 and to cor- 
rect them at higher temperatures, we undertook Monte 
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Carlo integration computations following accepted proce- 
dures [44j . However, the general sample- mean method, at 
first yields results with errors significantly too large at the 
small values of T* needed for ionic criticality. The rea- 
son is simply that the integrand of K3 3 is sharply peaked 
around the ground state, the peak sharpening as T* is 
lowered and hence becoming less frequently sampled. To 
improve the accuracy, we used a 'weighted sample-mean 
method', in which random numbers are generated with 
a weighting chosen to sample the integrand more often 
near the peak. Thus, in one dimension, for example, to 
evaluate I=J a f(x)dx, one needs to calculate 



We generalized this procedure to the geometry of the 
tetramer. For the radial integrals, we used the density 
function 



p(r) = A r exp(— X r r) with \ = C^^/T* 



(B3) 



and A r = 1/ J a exp(— A r r) dr. This weighting mimics the 
peaks in the integrand almost exactly. For the angular 
variables the optimal weighting is more complicated be- 
cause the peaks are Gaussian leading to an equation l|B2|l 
that cannot be simply inverted algebraically. Instead, we 
used exponential weighting 



/(Si) 
P{Xi) ' 



(Bl) 



where p(x) is the probability density function used to 
generate the random numbers, normalized in the interval 
[a, b], and n is the total number of random points xi S 
[a, b\. The weighted random numbers are generated from 
the uniform random numbers a £ [0, 1] by solving for x 
in 



P(x) = / p(x') dx' =a. 



(B2) 



The density function should be chosen so that this rela- 
tion can be solved for x algebraically. 



p(oj) = exp(-A w |w|) , 



(B4) 



with Lu — (p, 62 — 0i2—2tt/3 or #3 = 6*13— 2n/3, and with the 
normalizing integrals A~ x = J* exp(— A v |^|) and Ag 1 — 

f^ /3 exp(-A0|0|). We also chose X v = 1/(24V3T*) 1 / 2 

and A = 2.5/(24 V / 3T*) 1 / 2 so that the width of the peak 
in l|B4ll matched the width of the peak of the integrand. 
Finally, as a generalization of the Bjerrum procedure, we 
used a radial cut-off R = 0.196 a/T* which satisfactorily 
located the minimum of dK^ ^/dR. 

The results, which are well fit by (13 . X^tf) with the coef- 
ficients listed in Tabled agree closely with the consensus 
of the seventh order Pade approximants up to T* ~ 0.03; 
but they deviate strongly above T*=0.06: see Fig. 
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